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llJNif-lJrj: 


The  most  notable  changes  occurred  when  testing  the  ADI  program  on  the 
conditions  listed  on  page  100.   The  output  from  this  test  is  attached, 
and  it  shows  the  corrected  values  for  the  L2  norm  ratio  only.   The  L2  norm 
ratio  was  chosen  since  it  gives  the  test  indication  of  the  method's 
convergence  to  the  true  solution.   Although  tests  were  not  performed, 
discrepancies  in  the  results  appearing  on  pp.  102,  108,  and  110  would  be 
comparable  to  those  between  the  attached  output  and  page  100  of  the  report. 
The  results  of  the  model  problems  (pp.  96,  98,  10U,  106) ,  are  only  slightly 
affected  by  the  indexing  error  cited  above.   The  reason  being  that  in  the 
model  problems  described  in  the  report  the  region  is  homogeneous,  i.e.,  the 
coefficients,  Al,  are  a  constant  multiple  of  the  coefficients,  A2,  with  Al 
and  A2  being  constant  throughout  the  region. 


On  page  15,  in  equations  number  35  5  PE  and  PW  are  in  error.   They  should  be 
changed  as  follows : 


PE  =  -[ai(i,j)  +  a1(i+l,j)] 
PW  =  -[a^i.j)  +  ELjU-l.j)] 


r:bjm 
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1.   In  subroutine  HH1  (pp.  65,  66  and  again  on  pp.  86,  87)  the  index  Jl  was 
found  to  be  in  error.   The  remainder  of  the  program  following  INM1  =  N-l 
should  be  changed  as  follows: 


*       LI  =  1 

DO   5   I  =  1,N 

PW  =  Al(l+M)  +  A1(I+M+1) 

DO      k     J   =   II,INM1 
J   IS   THE  H   INDEX   fflD  Jl  IS   THE  Al   INDEX 

Jl  =  J+Ll 

PE   =  A1(J1)    +   A1(J1+1) 

PO   =   -(PE+PW) 

H(J)    =   -PE/(P0+W+PW*H(J-1)) 

PW  =  PE 
k      CONTINUE 

II  =  II+N 

INM1  =  INM1+N 

M  =  M+Nl 

LI  =  Ll+2 
?   CONTINUE 

RETURN 

END 


Note  1.   Where  an  *  precedes  a  statement,  a  statement  has  been  added  to 
the  program. 

Note  .2.   Where  **  precedes  a  statement,  a  statement  has  been  altered. 

As  one  can  see,  before  this  change  the  index  Jl  was  taking  on  the  values 
3  through  30,  then  33  through  60,  then  63  through  90,  etc...  .   However  Jl 
should  have  been,  first  3  through  30,  then  35  through  62,  then  6?  through  9^ 
etc. . .  . 

This  previous  change  in  the  ADI  program  does,  in  effect,  change  some  of  the 
results  noted  on  the  following  pages  of  the  report:  96,  98,  100,  102,  10 J), 
106,  108,  and  110. 
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1.     INTRODUCTION 


1.  1     The  Problem 


Let  Au  =  s  be  a  linear  system  of  algebraic  equations  resulting 
from  the  discretization  of  an  elliptic  boundary  value  problem.     A  di- 
rect method  for  computing  the  solution  is  the  same  as  factoring  A  into 
the  product    of  a  lower  triangular  matrix,    L,    and  an  upper  triangular 
matrix  U,    such  that  A  =  LU.      Because  L  and  U  are  not  sparse,    such 
methods  are  computationally  inefficient.     More  efficient  methods  are 
iterative,    not  direct.     A  particularly  efficient  iterative  procedure  is 
Alternating  Direction  Implicit  (ADI)  proposed  by  Peaceman  and  Rach- 
ford  [l],    and  Douglas  and  Rachford  [2].     See  Varga  [3]  and  Wach- 
spress  [4].      The  idea  of  this  method  is  to  partition  A  as  the  sum  of 
easily  invertible  matrices  H  and  V, 

(1)  A  =  H  +  V, 

and  to  compute  the  solution,    u,    as  the  limit  of  the  sequence  defined  by 

(2)  (H  +  wI)u(t+l/2)       =     S  -   (V  -  col)  u(t) 

(V  +  coI)u(t+1)  =     S-(H-coI)u<t+1/2) 

where  u  is  a  parameter  used  to  accelerate  convergence. 

Recent  iterative  procedures  proposed  by  Stone   [5]  and  Dupont, 

Kendall,    and  Rachford  [6],    are  based  on  finding  a  matrix  A  +  B      which 

a 


2 

can  be  factored  into  the  product  of  two  sparse  matrices,    L  and  U, 


(3)     A  +  B     =  LU. 


The  solution  of  (A+B  )  u  =  s  is  computed  as  the  limit  of  the  sequence 
defined  by 

(4)     (A  +  B    )  un+1   =  (A  +  B     )  un  -  T(Aun  -  s), 

where  a  and  j   are  parameters  to  accelerate  convergence.      The  pro^ 
cedures  of  Stone,    and  Dupont,    Kendall,    and  Rachford  differ  only  in 
the  method  used  to  determine  A  +  B   . 

a 


1.  2    Purpo 


se 


Dupont,    Kendall,    and  Rachford  succeeded  in  analyzing  their 
procedure,    obtaining  conditions  for  convergence  and  estimates  on  the 
rate  of  convergence.      Dupont  generalized  and  extended  this  analysis 
in  [7].     No  similar  published  information  is  available  concerning  the 

tone  procedure.      The  properties  of  this  procedure  have  been  derived 
principally  from  numerical  studies  as  described  by  Stone  in  [5],    al- 
though Stone  does  give  an  interesting  heuristic  analysis  in  the  same 
paper. 

The  purpose  of  this  thesis  is  to  study  the  Stone  procedure 
lly.      Solutions  of  the  Dirichlet  and  Neumann  boundary  value 
me   for   several  differential  operators  are  computed  by  an 


algorithm  based  on  the  Stone  procedure.     Solutions  to  the  same  prob- 
lems are  computed  using  ADI,    and  numerical  features  of  both  pro- 
cedures are  compared.      These  studies  provide  an  independent  veri- 
fication of  the  results  Stone  reported  in  [5], 


2.       MATHEMATICAL  PRELIMINARIES 

The  Dirichlet  problem  is  to  compute  the  solution  of 

(5)     --*-    (a    (x)    |H-).-^-t   (x)    &-)   =S(x),     x€D 

u(x)  =  f(x),   x  e  BD, 
where  x  =  (x   ,    x    ),    D  is  the  interior  of  a  compact  region,    and 
the  differential  operator, 

is  elliptic   in  D. 

The  Neumann  problem  is  to  compute  the  solution  of 


dx1  v  i      dx^    ax2  \  2     0x2  / 


=  S(x),       x€D 


O  x 


with  x,    D,    and  the  differential  operator  as  above. 

Let  D  be  the  interior  of  a  square  with  sides  parallel  to  the 
coordinate  axes  and  cover  D  with  a  square  grid  system  formed  by  N 
horizontal  and  N  vertical  lines,    replace  the  differential  operator  (5) 

i  ri  difference  operator,    and  replace  u(x)  with  a  vector  whose  conv 

2  2 

}  correspond  to  the  N     grid  points.      This  leads  to  a  set  of  N 

OUI   linear  equations  in  N     unknowns  where  each  equation 

uniquely  to  a  grid  point  of  the  system. 


Specifically,    let 


(7)     D  =  {(xrx2):     0  <  Xj  x2  <  lj.. 


Let  the  spacing  between  grid  points  be  uniform  and  defined  by 


(8)  h=  "nTT  ■  N>0' 


and  let 


(9)     xH  =  iA  Xj 

x2.   =  jA  x 


1  <  i  <  N        Ax     =  l/N 
1  <  j  <  N        Ax     =  l/N. 


Define  the  grid  system  D     over  the  region  D  by 

h 


(10)  D  =  {(ih,    jh)}. 

Denote  the  grid  point  (i  h,    j  h)  by  (i,    j).      The  boundary  dD     of  D,    is 
defined  by 


(ID  8Dh  = 


(i, 

0) 

0  <  i  <_N+1, 

(1> 

j) 

0  <  j  <  N+l, 

(i, 

1) 

0  <_j  <  N+l, 

(0, 

j) 

0  <  j  <  N+l. 

The  grid  system  is   shown  in  Figure   1.      Order  the  grid  points  of  D 

in  a  left-to-right,    down-to-up  fashion,    and  let  u  be  a  vector  whose 

components  u.    .  are  ordered  in  the  same  sequence  as  the  grid  points 
1»  J 


(0,1) 
(0,Nh) 


(0,h) 
(0,0) 


i 

(i,D 

r~ 

1 

i 
1 

i 
1 

i 

1 

1 

! 

S 

_ 

1       I 

1       1 
1       1       1 

1 

I 
I 

1 

1 
1 

1 

(h,0) 


(Nh.O)    (1,0) 


#   i 


Figure   1.      Grid  Point  System 


(i,j).   Thus,  u.  .  precedes  u.,  .  if  i  <  i',  and  u.  .  precedes  u    if 
1,  j  i  ,  J  i»  J  ^»  -t/ 

j  <  L 

Several  discretizations  of  the  expression 

are  possible.     For  example, 

(13)    h"2   {[a^i.j)  aj   (i+1,    j)]1/2  [u(i,  j)  -  u(i+l,  j)]  + 
[a^i,  j)  ax   (i-1,   j)]l/2  [u(i,  j)  -  u(i-l,  j)]  + 
[a2(i,  j)  a2  (i,  j  +  l)]l/2  [u(i,j)  -  u(i,j+l)]  + 
[a  ,(i,j)a,  (i,j-l)]l/2  [u(i,j)  -  u(i,j-l)]}, 


or 


h-2 
(14)  —    [a^i.j)  +  a^i+l.j)]  [u(i.j)  -  u(i+l,j)]  + 


[a^i,  j)  +  a^i-1,  j)]  [u(i,j)  -  u(i-l.j)]  + 


[a2(i,j)  +  a2(i,j+l)]  [u(i,j)  -  u(i,  j+1)]  + 


[a2(i,j)  +  a2(i,j-l)]  [u(i,j)  -  u(i,j-l)]. 


No  attempt  has  been  made  to  consider  the  mathematical 
merits  of  one  discretization  over  another.     The  criteria  used  result 
from  storage  problems,   not  mathematical  considerations.      The  need 


8 
to  minimize  the  storage  requirements  of  the  routines  used  in  this 
thesis  for  testing  the  ADI  and  Stone  algorithms  force  the  computa- 
tion of  the  elements  of  A,    at  almost  every  step  of  execution,    from 

the  coefficients,    a    (x)  and  a    (x).     Computing  the  elements  of  A 

2 
this  way  allows  the  storage  of  only  2N     +  4N  values  of  the  coeffi- 

2 
cients,    a^x)  and  a    (x),    instead  of  the  (at  most)  5N     -  2N  non-zero 

elements  of  A.      For  these  reasons,    the  discretization,    (14),    is  the 
most  practical  for  converting  the  boundary  value  problem  to  matrix 
form. 

The  exact  definition  of  these  elements  of  A  as  obtained  from 
(14)  depends  on  the  type  of  boundary  value  problem,    Dirichlet  or  Neu- 
mann.     Either  problem  may  be  written  in  the  form  Au  =  S,    where  A 

2  2  2  2 

is  an  N     x  N     matrix  if  Dirichlet,    and  (N+2)     x  (N+2)     if  Neumann. 

The  difference  in  the  order  of  A,    of  course,    is  due  to  the  fact  that 
for  the  Dirichlet  problem,    components  of  u  are  given  on  the  boundary, 
whereas  for  the  Neumann  problem,    they  are  not. 

If  (i,  j)  is  a  point  of  D     for  which  (i+_l,  j)  or  (i,  j+_l)  is  not  on 
the  boundary,    then  the  discretization  (14)  yields  the  following  values 
for  the  elements  of  A  for  either  the  Neumann  problem  or  the  Dirichlet 
problem: 


B.         =     -a_(ih,  jh)  -  a    (ih,  (j-l)h), 
i.J  2  2 

D.    .     =     -a.(ih,  jh)  -  a   (  (i-l)h,   jh), 
i.J  1  l 

(15)      E.    .     =     -(B.    .  +  D.    .  +  F.    .  +  H.    .), 

F.         =     -a.  (in,  jh)  -  a    (  (i+l)h,   jh), 
i.J  1  ! 

H.    .     =     -a_(ih,  jh)  -  a_(ih,    (j  +  l)h), 
i.J  2  2 

and  S.    .  =  f(i,  j). 
i»  J 

The  treatment  of  what  happens  at  the  boundary  is  straight- 
forward.    Let  (i,  j)  be  a  grid  point  for  which  an  adjacent  point  is  on 

the  boundary  of  D,  ,    that  is,    i  =  1  or  N  or  j   =  1   or  N.     Suppose  i  =  1, 
h 

1  <  j  <  N.      The  other  cases  are  the  same.      For  the  Dirichlet  problem, 


(16)    D     .  =  0 
i.J 


and 


(17)  S.    .  =  f(i,  j)  +  [a    (i,  j)  +  a    (i-l),j)]  u(i-l,  j). 

The  other  elements  are  defined  as  in  (15).      This  completes  the  de- 
scription of  the  Dirichlet  problem.     For  the  Neumann  problem  the 

elements  of  A  and  the  inhomogeneous  part,    S       ,    are  defined  exactly 

i.  J 

as  in  (15). 

Now,    suppose  i  =  0,    1  <  j  <_N.      Only  the  discretization  cor- 
responding to  the  Neumann  problem  is  defined.     We  have, 


(18)    D.    .  =  0 


and 


(19)  f(i,j)  =  -2  [a^i.j)  +  a^i+l.j)]. 
The  other  elements  are  defined  as  in  (15). 
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3.       STONE  ALGORITHM 
3.  1     The  Factorization 

The  factorization  of  A  into  the  product  of  a  lower  triangular 
matrix,    L,  by  an  upper  triangular  matrix,    U,    such  that  A  =  LU,   does 
not  yield  a  sparse  fadtorization.     Stone  [57]  suggests  constructing  L 
and  U  to  be  sparse  matrices  with  three  non-zero  diagonals,    each  of 
which  coincide  with  appropriate  diagonals  of  A.      The  product  of  L 
and  U  now  defines  an  altered  form  of  A, 


(20)    A  +  B     =  LU, 
Of 

which  makes  the  direct  solution  of  (A  +  B    )  u  =  s  computationally  easy. 


Let  the  elements  of  L  and  U  be  defined  by 


(21)  (Lu).    .  =  b.    .  u.    .  +  c.    .  u.    ,     .  +  d.    .  u.    ., 
i,j         i,j     i,j-l  i,j     i+l, j         i,  j     i,j 

(Uu).    .  =  u.    .  +  e.    .  u.    ,     .  +  f.    .  u.    .    ,. 

1>i         !>J         i.J  i+l, J         i,  J     i.j  +  1 


The  product  LU  is  defined  by 


(22)  [(LU)u].    .  =  b.    .  u.    .    ,   +  b.    .  e      .  u.    , 

i,j         i,j      i,  j-1  l,  j     i,j     i+l.j-l 

+  c        u  .  +  (d.    .  +  b.    .  f.    .       +  c.    .  e.    .     .)  u.    . 

X»J     1_1'J  i,J  i*J     i»  J-1  i,J     i-l,  J       1,3 

+  d.    .  e.    .  u.         .  +  c.    .  f.    ,    .  u;    , 

1,3     1,3     1+1,3         i,J    i-l, J     i-l,3+! 

+  d.    .  f.    .  u.    .    ,. 

1,3    1,3     1,3  +  1 

Stone  computes  the  elements  of  L  and  U  by  the  following  al- 
gorithm: 
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b.    .  =  D.    .  -  ry  G.    ., 
i.J  LJ  i.J 


c.    .   =  D.    .  -  a  C.    ., 
l»J  LJ  i.J 


(23)     d.    .  =  E.    .+a  (C.    .  +  G.    .)  -  b.    .  f.    .   .   -  c.    .  e.    ,    ., 

LJ  LJ  i.J         i.J  i.J    i.J-1         i.J     i-l.  J 


»  J 

F. 
i 

.,)' 

•Of 

c. 

1, 

j 

e . 
i, 

d. 

-.j 

. 

f. 
if 

J 

H. 
h 

j" 

a 

G. 

if 

j 

d. 

9 

I.J 

where        0  <   q,  <  1. 

3.  2    The  Iterative  Method 


The  iteration  defined  in  (4), 


(24)  (A  +  N)un+1   =  (A  +  N)un  -  T   (Aun-s), 


is  computed  as  follows  (Stone  [5])  for  j   =  1.      Define 


(25)    Rn  =  s   -  Au11 


and 


/_/  .    pn+l  n+1  n 

(26)  6  =  n  -  u   . 


Solve 


(27)  LV  =  Rn 


for  V  and 
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(28)  U  6n+1  =  V 

n+1        .  , 
for  5        .     Then  set 

._.,     n+1  n        n+1 

(29)  u  =  u     +  6 

3.  3   Application  Notes 

Stone  [5]  suggests  a  variation  for  alternate  steps  of  the  itera- 
tion which  has  been  followed  in  this  thesis.      For  all  odd  iterations  be- 
ginning with  the  first,    the  solution  is  obtained  at  the  (i,  j)  grid  points 
by  starting  with  the  j  =  0  row  and  i  taking  on  the  values  0  <  i  <_  N. 
Then  j  is  incremented  by  1  and  i  again  assumes  the  above  range  of 
values,    until  j  =  N.      For  all  even  iterations  beginning  with  the  second, 
this  process  is  reversed,    that  is,    the  solution  is  obtained  at  the  points 
in  the  j   =  N  row  with  i  assuming  the  values  0  <  i  <_  N.      Then  j  is  de- 
cremented by  1  and  i  again  assumes  the  range  of  values  0  <  i  <  N. 

At  the  even  iteration  steps  then,    the  computation  of  the  solu- 
tion begins  at  the  upper  left-hand  grid  point,    (1,  Nh),    and  moves  in  a 
left-to-right,    up-to-down  order  over  the  grid  instead  of  the  left-to-right, 
down-to-up  order  of  odd  iterations.      This  new  ordering  of  the  grid  points 
necessitates  recomputing  the  coefficients  of  the  altered  matrix,    A  +  N, 
for  the  even  steps. 

A  detailed  description  of  the  Stone  program  is  provided  in  Ap- 
pendix A. 
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4.       API  ALGORITHM 

Let  H  and  V  be  discretizations  of 


(30,-&:(ai(x)£:) 


and 


<31>-i^(a2<*>iij)- 

respectively  [8],    so  that 

(32)  A  =H  +  V. 

The  solution  u  of  Au  =  S  is  computed  as  the  limit  of  the  sequence  de 
fined  by 

(33)  (H  +  co     ,    I)un+1/2  =  S  -  (V  -  co         I)un, 

n,  h  n, v 

(V  +  co  I)un+1        =  S  -  (H  -  co     u  I)un+1/2. 

n, v  n, h 

Since  H  and  V  are  real,    tridiagonal,    symmetric  matrices,    both 

(H  +  co  I)  and  (V  +  co  I)  are  positive  definite  and  their  inverses 

n,  h  n,  v 

exist.      The  ADI  algorithm  is  based  on  (33). 

A  detailed  description  of  the  ADI  program  is  provided  in 
Appendix  B. 
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5.       CODING  OBJECTIVES 

5.  1    Storage  Minimization 

Let  A  be  the  matrix  formed  by  discretizing  the  differential 
operator 

(34)    -&_/a    (X)   rM+    -*-fa,(x)   -M. 

bx:  \  i      dXl  ;     ax2  \  z      ax2 ; 

Let  (i,  j)  be  a  grid  point,    and  let 

PS  =  -[a2  (i,  j-1)  +  a2  (i,j)], 

PE  =  -[a     (i-l,j)  +  ax   (i,  j)], 

(35)  PW  =  -[ax  (i,j)  +a2  (i+l,j)], 
PN  =  -[a2  (i,j)  +  a2  (i,j+l)], 
PO  =  -(PS  +  PE  +  PW  +  PN). 

If  the  discretization  in  (14)  is  followed,    the  row  of  A  corresponding 
to  the  grid  point,    (i,  j),    1  <   i,j  <  N,    is 

(36)  PS  .    .    .    PE      PO      PW  .    .    .    PN, 

2  2 

where,    if  the  order  of  A  is  N     x  N    ,    PS  and  PW  are  N  positions  re- 
moved from  PO  on  the  main  diagonal. 

All  possible  values  of  PS,    PE,    PO,    PW,    and  PN  form  five 

diagonals.      The  maximum  number  of  non-zero  elements  of  these 

,.  .      ,  2  2         2         2  2 

diagonals  are  respectively,    N     -  N,    N   ,    N   ,    N     and  N     -  N,    totaling 
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5N     -  2N  non-zero  elements  of  A.     Actually,   values  of  PE  corres- 


ponding to  i  =  1  are  zero.     Accordingly,    the  number  of  non-zero  ele- 

2 
ments  of  the  diagonal  formed  by  the  possible  values  of  PE  is  N     -  N. 

Since  A  is  symmetric,    the  number  of  non-zero  elements  in  the  diagonal 

2  2 

formed  by  PW  is  also  N     -  N.      Therefore,    the  bound  5N     -  2N  may  be 

replaced  by  5N     -  4N. 

As  a  result,    storage  of  A  is  possible  with  only  5N     -  4N  loca- 
tions,   although  the  treatment  of  the  boundary  points  required  to  obtain 
this  bound  is  a  somewhat  hazardous  programming  feat. 

By  using  the  fact  that  the  diagonals  of  A  are  computed  from 
the  coefficients,    a    (x)  and  a    (x),    it  is  possible  to  reduce  the  storage 
requirements  even  further,    although  at  the  expense  of  program  ef- 
ficiency,   by  storage  of  the  necessary  values  of  a    (x)  and  a    (x).     The 
necessary  values  of  a    (x)  are  values  a    (i,  j)  for  which  0  <_i  <  N  +  1, 
1  <_  j  <_  N,    and  the  necessary  values  of  a    (x)  are  all  values  for  which 

1  <  i  <N,    0<j^N  +  l   (see  Figure  2).     A  total  of  2JN  (N  +  1)]  = 

.     2 

2N     +  2N  storage  locations  are  therefore  required  to  store  the  coef- 

2 
ficients,    a    (x)  and  a    (x),    a  reduction  from  the  5N     -  4N  locations  re- 
quired to  store  A. 

The  advantage  of  coding  the  Stone  and  ADI  procedures  using 
only  the  coefficients  a    (x)  and  a    (x)  is  simply  that  the  size  of  the 
largest  problem  which  fits  in  core  is  thereby  increased,    allowing 
asymptotic  studies  of  the   rate  of  convergence.      The  disadvantage,    of 


17 


r 


(0, 

1) 

• 

(0, 

Nh)p 

1 

1- 
1 

1 

1 

r 

i 

i 
i 

i 

i- 

I 1 T 

I  I 


s.(0,h)    L_ 


i — T""~n 


(i,D 


(0,0) 


i I L__L 

(h,  0) 


JNh.Nh) 

I 


1 

I 

--I 
I 
I 

1 

I 

-H 

I 

I 
i 

I 


_J I 

(Nh,0)    (1,0) 


L 


V 

ax(x) 

igure  2.       Definition  of  Differential  Equation  Coefficients,    a    (x)  and 
a    (x),    and  Their  Relation  to  the  Grid  Point  System 


18 
course,  is  that  computing  PS,  PE,  etc.  ,  at  each  step  of  the  algorithm 
grossly  reduces  efficiency. 

5.  2    Reducing  Execution  Time 

An  effort  was  made  to  localize  all  multiplications  and  divisions 
into  successive  additions  and  subtractions  as  Gear  [9]  suggested  in  his 
paper  on  compiler  techniques.      For  example,    the  calculations 

DO  30        J=l,N-2 

PE     =     Al(J*N+4)  +  Al(J*N+5) 
PN    =    A2(J*N+1)  +  A2((J+J)*N+1) 
PW    =     Al(J*N+3)  +  Al(J*N+4) 
(37)     PS     =     A2((J-1)*N+1)  +  A2(J*N+1) 
PO    =     -  (PE+PN+PW+PS) 

RN(J*N+1)  =  PS*T((J-l)*N  +  l)+PO*T(J+N+l) 
+  PE*T(J*N+2)+PN*((J+l)+N+l) 


in  subroutine  PRODMT  are  localized  as  follows: 


19 


Ml  =  N+l 

M2  =  N 

DO  30         J=N+1,   N**2-2*N+1,   N 

Ml  =  Ml+3 

M2  =  M2+1 

(38)      PE  =  A1(M1)+A1(M1  +  1) 

PN  =  A2(M2)+A2(M2+N) 

PW  =  A1(M1-1)+A1(M1) 

PS  =  A2(M2-N)+A2(M2) 

PO  =  -(PE+PN+PW+PS) 

RN(J)  =  PS*T(J-N)+PW*T(J)+PE*T(J+1)+PN*T(J+N) 


Execution  time  has  been  saved  by  replacing  the  multiplication  J*N  in 
(37)  by  successive  additions  in  (38). 
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6.       NUMERICAL  RESULTS 

6.  1    Description  of  Problems 

Four  different  problems  are  considered  in  this  thesis.      Each 

was  solved  over  a  square  grid  system  with  30  points  on  a  side.     Some 

results  were  obtained  for  10x10  and  20x20  grid  systems,   but  these  are 

not  reported  here  since  they  were  essentially  the  same  as  those  for  the 

30x30  grid.      The  solution,   u(x),    of  (5)  and  (6)  was  arbitrarily  chosen  to 

be  everywhere  equal  to  unity.     The  right-hand  sides,    S.    .,   were  com- 

1»  J 

puted  from  (Au).    .  =  S.    .,   and  then  used  to  solve  (Au).    .  =  S.    .  for  u.    .. 
i.J         if  J  l,  j         l,  j  l,  j 

All  initial  vectors  were  everywhere  equal  to  zero,    except  for  compo- 
nents  1,    2,    N  +  1,    and  N  +  2  which  were   -1,    +1,    +1,    and  -1  respec- 
tively. 

The  first  problem  considered  is  the  model  problem,   with  a    (x) 
and  a    (x)  both  equal  to  -.  5  over  the  entire  region.      The  coefficients  of 
A  were  obtained  by  the  discretization  (14)  of  the  operator  (34). 

The  second  problem  considered  is  the  generalized  model 
problem  in  which  the  ratio  of  a    (x)  to  a    (x)  was  equal  to  100  with 
a    (x)  =  -.  5  and  a    (x)  =  -.  005. 

The  third  problem  was  taken  from  the  problem  described  by 
Stone  [5].  Refer  to  Figure  3.  In  region  A,  the  ratio  of  a  (x)  to  a  (x) 
is  equal  to  one.      This  ratio  in  region  B  is  equal  to  ten.      In  region  C 

ratio  is  one -tenth,    and  in  region  D,    one.     However,    in  D  the  values 
of  'i    (x)  and  a    (x)  are  different  by  an  order  of  ten  from  those  in  region 
A. 
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Figure  3.       Heterogeneous  Test  Problem  Geometry 
(Region  Defined  on  Unit  Square) 
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The  last  problem  is  the  same  as  the  third  with  one  change. 
In  region  A,  a  (x)  and  a  (x)  were  varied  in  a  random  manner  by  the 
aid  of  a  random  number  generator.     The  range  of  variation  satisfied 

(39)  0  <  a^x),    a2(x)  <   1. 

These  four  cases  are  taken  from  those  presented  by  Stone  [5], 
and,    in  this  sense,    the  results  reported  herein  independently  verify  his 

studies. 

6.  2    Iteration  Parameters 
6.  2.  1     Stone  Program 

The  suggestions  of  Stone  [5]  were  followed  in  the  selection 
of  parameters  for  the  computations.     A  set  of  nine  parameters  was 
used  in  each  of  the  four  test  cases  for  both  the  Dirichlet  and  Neumann 
problems.      Each  parameter  was  used  twice  in  succession,    yielding  a 
total  of  18  parameters  per  cycle.      This  set  of  18  parameters  was 
broken  into  three  sets  of  6  each,    where  each  set  of  6  covered  the 
range  0  <_    #  <  1 .      For  example,    if  a     <  0Ly  <   •  •  •    <   aQ  are  the 
parameters,    then  they  are  applied  to  the  iterative  method  in  the 
following  sequence:    a    ,    oy    «6»    %>    0>y    Qy    Of g>    ag>    OL^,    a 5>    Oi 2> 

"/,   ®7'   °V   a4-    &4,   ofj,   cvj. 

The  minimum  parameter  was  always  zero  and  the  maxi- 
mum parameter  was  computed  by 


(40)   a  =  1-min 

max 
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2  A  x 


1 


2  A  x. 


1  + 


a    (x)  A  x 


2* 


L     1  + 


a   (x)  A  x. 


a^x)  A  x   ' 


*2(XMX1 


The  intermediate  parameters,    0  <    <y      </y  ,   m  =  l,    ...,M-1, 

"~     in  ■■     max 

were  computed  from 


(41)  ^  =  1  -  (1  -  ttmax) 


m 
M-l 


,       m  =  l,    ...,    M-l, 


where  M  =  9.     The  value  of  M  is  the  number  of  parameters  in  a  cycle 
of  18  iterations. 


6.  2.  2    ADI  Proi 


ram 


Let  u       be  the  approximation  to  the  solution  of  Au  =  S 


computed  by  the  ADI  procedure  defined  in  (2).      In  terms  of  u 


2n-2 


(42)  u        =  T     u 

c.n  n     cx\.~  £. 


where 


(43)  T      (V  +  co  I)"1   (H  -  co  I)  (H  +  co  I)"1  (V  -  o>  I), 

n  v,  n  h,  n  h,  n  v,  n 


Accordingly, 


n 


<«>  ll»2n  h  <_  II  TT  t  ||2  ||»0n2. 

i=0 


Assume  that  V  and  H  commute  and  that  they  have  the  same  eigenvectors 
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Let  x  be  such  an  eigenvector  with  eigenvalues   \  and  ^  corresponding 
to  H  and  V  respectively.     Then  x  is  an  eigenvector  of 


n 


77  t. 


1=0 


with  corresponding  eigenvalue 


n 


(45) 


TT 


a  -  w  X  -  w,     • 

**         v,  1  h,  1 


1=0    **         v,  1  h,  1 


Conversely,    every  eigenvalue  of 


n 


IT  Ti 


1=0 


is  of  the  form  (45)  where  ^  and  \   are  eigenvalues  of  H  and  V. 

Ideally,    the  parameters  ax     .  and  go      .  should  be  chosen 
7  r  h,  1  v,  1 

to  minimize  the  absolute  value  of  the  expression  in  (45)  when  ^  and  X 
range  through  all  the  eigenvalues  of  H  and  V  respectively.     Such  a 
choice,    of  course,    minimizes 

n 


TT  T; 


1=0 


since  this  quantity  is  equal  to  the  largest  eigenvalue.      Parameters,    how- 
ever,   are  selected  somewhat  differently.      Let  a  and  b,    c  and  d,    be 
bounds  on  the  eigenvalues  of  H  and  V,    i.    e.  ,    0  <  a  <  \  <  b  and 
0  <c   <   n  <   d.      The  bounds   a  and  b  are  not  positive  for  the  Neumann 
>blem.      For  the  testing  performed  in  this  thesis,    parameters  a)      . 


and  co      .  were  chosen  to  minimize  the  rational  expression 
v,  1 
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max 
(46)*<W<b 
c  <  Z  <  d 


n 


TT 


W  -  co,     .  Z    -  co      . 

h, l  v,  l 


i=o     W  + 


Vi 


Z  +  co 


v,  l 


according  to  an  algorithm  of  Wachspress  [4]. 

For  the  Neumann  problem,    a  =  c  =  0  and  for  any  choice 
of  parameters,    the  expression  in  (46)  is  equal  to  1.      Therefore,"  opti- 
mal parameters  for  the  Neumann  problem  must  be  chosen  more  care- 
fully.    In  order  to  make  this  choice,    assume  that  the  null  spaces  of  H 
and  V  coincide  with  the  null  space  of  A,    A  resulting  from  the  Neumann 
problem.      The  null  space  of  A  is  the  subspace  spanned  by  any  non-zero 
constant  vector.      Let  N  denote  the  null  space  of  A  and  N     its  orthogonal 
complement.      The  restriction  of 

n 


TTTi 


1  =  o 


to  N  is  singular  whereas  the  restriction  of 


n 


U  Ti 


1  =  o 


to  N     is  not.      The  eigenvalues  of 


n 


TT  Ti 


i  =  o 


restricted  to  N      have  the  form  of  the  expression  in  (45),   where  the 


26 
lower  bounds  on  the  eigenvalues  are  positive. 

Let  a  and  c  be  the  lower  bounds.      The  smallest  eigen- 
values of  H  and  V  are,    of  course,    zero.     Suitable  values  for  a  and  c, 

for 

n         _ 

TT    Ti 

1=0 

restricted  to  N   ,    are  the  next  smallest  eigenvalues  of  H  and  V.     Since 

a  and  c  are  positive,    a  minimizing  sequence  of  parameters  co         and  00 

h,  l  v,  1 

reduces  the  norm  of 

TT     Ti 

i  =  o 
when  restricted  to  N   . 

The  method  followed  here  to  select  parameters  for  the 
Neumann  problem,  therefore,  is  to  determine  the  next  smallest  eigen- 
values of  H  and  V.      Then,    using  these  as  lower  bounds  in  (46),    select 

a  sequence  to  minimize  (46).      Let  e     be  an  initial  error  vector  and 

o 

write 


(47)  e      =  n  +  n 
o 

where  n  f  N  and  n      f  N    .      The  minimizing  sequence  therefore  reduces 

to 

n 


(48)    I     I       T.n   . 
i  =  o 
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(49)  T.  n  =n, 

1 
l  =  o 


the  error  e     is  reduced  to 
o 

n  n 


(50)"]  T.  e     =  n  +  T        T.  n"1, 

!    '       i    o  :    '       i 

1  =o  1  =  o 

which  is  approximately  a  vector  in  the  null  space  of  A. 

For  the  tests  reported  in  this  thesis,    a  power  method 

[10]  was  used  to  determine  the  smallest  eigenvalues  of  H  and  V  by  deter- 

-1  -1 

mining  the  largest  eigenvalues  of  H       and  V      .      For  the  method  to  work, 

successive  approximations  to  the  eigenvector  corresponding  to  the 
smallest  eigenvalue  must  not  lie  in  the  null  space  of  H  or  V.        This  con- 
dition holds  for  succeeding  approximations  if  the  first  approximation 
is  in  the  orthogonal  complement  of  the  null  space  of  H  and  V,    for  H 
and  V  are  symmetric  and  preserve  orthogonal  complements.      To  ob- 
tain an  initial  vector  orthogonal  to  the  null  space  of  H  and  V,    let 

u,    =  -1,    u     =  1,    u  =  1,    and  u_.    ,  =  -1.     All  other  components  are 

1  2  N+l  N+2  r 

zero.      It  is  simple  to  verify  that  this  vector  is  orthogonal  to  any  vector 
which  is  constant  on  horizontal  or  vertical  grid  lines,    and  is  therefore 
orthogonal  to  the  null  spaces  of  H  and  V  respectively. 
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6.  3    Comparison  of  Iterative  Methods 

Figures  4-7  and  12-15  correspond  to  each  of  the  four  test 
cases.      The  ordinate  is  the  absolute  value  of  the  maximum  residual 
error  in  (25),    normalized  by  dividing  by- the  sum  of  all  positive  right- 
hand  sides,    S.    ..      The  abscissa  is  the  computational  work  necessary 
to  achieve  the  given  maximum  normalized  residual.      One  unit  of  com- 
putational work  is  equivalent  to  one  iteration  of  the  Stone  procedure  or 
approximately  1.44  ADI  iterations.     Since  the  plotting  of  results  after 
each  iteration  causes  an  oscillating  curve  (due  to  parameter  cycles), 
only  the  relative  minimums  were  considered,    and  a  straight  line  was 
used  to  approximate  these  points  in  such  a  way  as  to  facilitate  visual 
comparisons  of  the  rates  of  convergence. 

Figures   8-11  and  16   -   19  also  correspond  to  each  of  the  four 
test  cases  with  both  boundary  conditions.      The  ordinate  in  each  figure 
is  the  L,     relative  error  defined  by 

ll«n  - » ||2 

(50) 


n  o,. 

u      -u    || 


The  abscissa  is  the  corresponding  computational  work.     Again  a  straight 
line  approximation  fitted  to  the  relative  minimum  points  of  the  graph 
was  used  to  simplify  the  figures.     See  Appendix  D  for  the  numerical 
data. 


< 

D 

a 

w 

W 

X 

Q 

W 

N 
i— i 

< 

a; 

o 

x 
< 


-3 

10 

io-4 

\ 

■           Stone 
■  ADI 

io-5 

\^\ 

-6 
10 

\       ^ 
\ 

.    -7 
10 

\ 

. 

io"8 

\ 

\ 

ln-9 

\ 

t                     1 

» 

29 


0 


10 


20 


30 


40 


COMPUTATIONAL  WORK 
(NO.    OF  ITERATIONS  FOR  STONE  PROCEDURE) 


Figure  4.     Model  Problem  -  Dirichlet  Boundary  Conditions 
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Figure  5.      Generalized  Model  Problem  -Dirichlet  Boundary  Conditions 
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Figure  6.     Heterogeneous  Problem  with  Homogeneous  Subregions   -  Dirichlet 
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Figure  8.     Model  Problem  -  Dirichlet  Boundary  Conditions 
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Figure  9.      Generalized  Model  Problem  -  Dirichlet  Boundary  Conditions 
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Figure  10. 
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Figure  12.     Model  Problem  -  Neumann  Boundary  Conditions 
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Figure  13.     Generalized  Model  Problem  -  Neumann  Boundary  Conditions 
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Figure   16.      Model  Problem  -  Neumann  Boundary  Conditions 
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Figure   17.      Generalized  Model  Problem  -  Neumann  Boundary  Conditions 
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6.  3.  1     Dirichlet  Problem 

Figures  4  and  8  present  data  for  the  model  problem  with 
Dirichlet  boundary  conditions.     Note  that  ADI  converges  at  a  faster  rate 
than  the  Stone  method. 

Data  for  the  generalized  model  problem  with  Dirichlet 
boundary  conditions  are  presented  in  Figures  5  and  9.     ADI  maintains 
approximately  the  same  rate  of  convergence  as  in  the  model  problem; 
however,    there  is  a  significant  increase  in  the  rate  of  the  Stone  method. 

Figures  6  and  10  exhibit  the  results  for  the  heterogeneous 
problem  with  homogeneous  subregions  with  Dirichlet  boundary  conditions, 
The  complexity  of  this  and  the  next  test  case  is  comparable  to  the  com- 
plexity of  practical  problems. 

The  last  case  for  the  Dirichlet  problem  is  the  heterogen- 
eous region  with  random  variations  in  subregions,    and  the  data  for  this 
problem  are  presented  in  Figures  7  and  11.     As  in  the  previous  case, 
the  Stone  method  converges  faster  than  ADI.     Also,    the  ADI  rate  of  con- 
vergence decreased  as  the  test  cases  increased  in  complexity,    whereas 
the  rate  of  convergence  of  the  Stone  method  has  either  increased  (gen- 
eralized model  problem)  or  remained  nearly  the  same. 

Table   1   indicates  the  work  required  by  either  method  to 
achieve  a  maximum  normalized  residual  of  10"°  for  the  four  test  cases. 
Where  required,    straight  line  extrapolation  was  used  to  obtain  the  data. 
Note  that  except  for  the  generalized  model  problem,    the  Stone  method 
was  less  sensitive  to  the  nature  of  the  problem  itself. 
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method 


Stone 


ADI 


TABLE  1 


Dirichlet  Problem 


Work  Requirements  of  Iterative  Methods 


model 


31 


13 


RATIO: 
ADI/STONE      0.  42 


generalized 
model 


1.  60 


heterogeneous  model 
homogen.  random  no. 

subregions  subregions 


27 


34 


1.  26 


28 


85 


3.  12 


6.  3.  2     Neumann  Problem 

Figures   12  and  16  present  the  data  for  the  model  problem 
with  Neumann  boundary  conditions.     Again  ADI  converges  at  a  faster 
rate  than  Stone.      Note,    however,    that  for  the  Neumann  conditions,    both 
rates  have  decreased  slightly. 

For  the  generalized  model  problem  with  Neumann  bound- 
ary conditions,    Figures   13  and  17  show  that  the  Stone  algorithm  is  con- 
verging faster  than  ADI.      For  either  method,    convergence  is  slower 
than  for  Dirichlet  boundary  conditions. 

Figures   14  and  18  present  the  data  for  the  heterogeneous 
problem  with  homogeneous  subregions.      In  this  case,    with  Neumann 
boundary  conditions,    the  faster  rate  of  convergence  is  again  obtained 
by  the  Stone  algorithm. 

I "inally,    the  heterogeneous  model  with  random  subregions 
esented   in   Figures   15  and   19  also  shows  the  superiority  of  the  Stone 

jorithm.     Convergence  rates  have  consistently  decreased  for  both 
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algorithms  with  the  increasing  complexity  of  the  test  cases. 

Table  2  is  a  summary  of  the  work  required  by  both 

-5 
algorithms  to  reach  a  maximum  normalized  residual  of  10       for  the 

test  cases  with  Neumann  boundary  conditions.      This  table  is  patterned 

after  a  table  of  Stone   [5],    and  straight  line  extrapolation  was  used  where 

necessary.      It  is  apparent  that  the  Stone  algorithm  is  less  sensitive  to 

the  nature  of  the  test  case. 


TABLE  2 


Neumann  Problem 


Work  Requirements  of  Iterative  Methods 


generalized 

method 

model 

model 

Stone 

57 

33 

AD  I 

27 

46 

RATIO. 

ADI/STONE 

0.48 

1.41 

heterogeneous  model 
homogen.  random  no. 

subregions  subregions 

41  38 

45  66 


1.  10 


1.  74 


6.  4  Summary 

(i)  The  Stone  algorithm  was  more  effective  than  ADI  as  the  com- 

plexity of  the  problems  being  solved  increased,    although  the 
parameters  used  for  ADI  were  optimal  and  no  method  is 
known  for  determining  optimal  Stone  parameters.     ADI  rates 
of  convergence  tended  to  decrease  with  increasing  problem 
complexity  while  Stone  rates  were  relatively  stable. 
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(ii)       Convergence  rates  for  both  methods  were  faster  for  problems 
with  Dirichlet  boundary  conditions  than  for  Neumann  boundary- 
conditions . 
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7.    CONCLUSION 

In  view  of  the  results  obtained  by  this  author,    the  Stone 
method  is  in  fact  superior  to  ADI  methods.     Most  importantly,    the  re- 
sults obtained  by  Stone  [5]  have  been  verified.     It  is  imperative  to  note, 
however,    that  since  the  problems  considered  in  this  paper  were  not  ex- 
actly the  same  as  Stone's,    the  numerical  data  did  fluctuate  somewhat. 
Nevertheless,    the  results  were  still  verified;  furthermore,    the  compari- 
sons for  both  Dirichlet  and  Neumann  boundary  conditions  further  solidifed 
them. 
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APPENDIX  A 
STONE  PROGRAM  (DIRICHLET) 
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C STONE  PROCEDURE 

C FOR  ITERATIVE  SOLUTION  OF  IMPLICIT  APPROXIMATIONS  OF 

C TWO  DIMENSIONAL  PARTIAL  DIFFERENTIAL  EQUATIONS 

C 

C FOR  THE  LINEAR  SYSTEM  OF  EOUATIONSt  A*T»Q,  THIS  PROGRAM 

C COMPUTES  THE  MATRIX  A+N  WHICH  IS  FACTORED  INTO  THE  PRODUCT 

C OF  TWO  SPARSE  MATRICES,  L  AND  U,  SUCH  THAT  A+N=L*U.  THE 

C SOLUTION  OF  (A+N)*T*Q  IS  COMPUTED  AS  THE  LIMIT  OF  THE  SEQUENCE 

C DEFINED  BY 

C ( A+N ) *T ( N+ 1 ) « ( A+N ) *T ( H )-TAU* ( A*T  <  N )-Q ) , 

C WHERE  TAU  IS  A  PARAMETER  TO  ENHANCE  CONVERGENCE. 

C 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  AK960),A2(960),Q(900),RN(900),T<900),BL(870),CL(900>, 
XDL(900),EU(900),FU(900) ,PARAM(9) 

C INPUT  M,A1,A2,Q,T,N,ITMAX 

C 

C NOTE:  IT  IS  NECESSARY  FOR  THE  USER  TO  INPUT  THE  FOLLOWING 

C VARIABLES  INTO  THIS  PROGRAM: 

C M  =  ONE  HALF  THE  NUMBER  OF  PARAMETERS  PER  CYCLE  (THERE 

C ARE  2*M  PARAMETERS  PER  CYCLE  CALCULATED  IN  SUBROUTINE 

C PARAMS). 

C A1,A2  =  THE  VECTORS  OF  DIFFERENTIAL  EQUATION  COEFFICIENTS. 

C 0  =  THE  VECTOR  CONTAINING  THE  RIGHT  HAND  SIDES  OF  THE 

C LINEAR  SYSTEM. 

C T  =  THE  INITIAL  VECTOR. 

C N  =  NUMBER  OF  GRIDPOINTS  ON  ONE  SIDE  OF  THE  SOLUTION 

C GRID.  (THERE  ARE  N**2  GRIDPOINTS  IN  THE  TOTAL  SQUARE 

C GRID  SYSTEM). 

C ITMAX  =  TERMINATION  CRITERIA  OR  THE  NUMBER  OF  ITERATIONS 

C TO  BE  PERFORMED. 

C 

C THE  INPUT  THAT  FOLLOWS  IS  A  SAMPLE  OF  THE  NECESSARY  DATA. 

ITMAX=32 

N  =  30 

M  =  9 

DO  109  1=1,960 

Al ( I )=-.5D00 

109  A2( I  )=-.5D00 
DO  110  1=1,900 
T( I  )=0.0D00 

110  0(1  )=0.0D00 

C CALCULATE  THE  PARAMETERS  FOR  A  CYCLE  OF  2*M  ITERATIONS  FOR 

C CONSTANT  OR  VARIABLE  Al  AND  A2. 

CALL  PARAMS(N,A1,A2,M,RN,NSQP2N) 

C it  IS  UP  TO  THE  USER  TO  ARRANGE  THE  PARAMETERS  AS  HE  WANTS  THEM 

C USED  IN  THE  CYCLE  IN  THIS  SECTION.  THE  PARAMETERS  HERE  ARE 

C ARRANGED  IN  THE  FOLLOWING  MANNER.  IF  THE  PARAMETERS  ARE 

r- RN(  1  )<RN(2X...<RN(9)  ,  THEN  THIS  ORDER  IS  EQUIVALENT  TO  THE 

C SEQUENCE  RN ( 9 ) , RN ( 9 ) , RN ( 6 ) , RN ( 6 ) ,RN ( 3 ) ,RN ( 3 ) , RN ( 8 ) , RN ( 8 ) ,RN ( 5 ) , 

C RN(5)  ,RN(2)  ,RN(2)  ,RN(7)  ,RN(7)  ,RN(4)  ,RNU),RN(  1),RN(  1)  . 

PARAM( 1 )=RN(9) 

PARAM(2 )=RN(6) 

PARAM( 3  )=RN(3) 

PARAM(4)=RN(8 ) 

PAfcAM( 5)=RN(5) 

PAPAM( h) =RN(2 ) 
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PARAM(7)=RN(7) 

PARAM(8)»RN(4) 

PARAM(9)»RN(1) 
—COMPUTE  THE  SUM  OF  THE  POSITIVE  RIGHT  HAND  SIDES  FOR  USAGE  IN 
—THE  NORMALIZING  OF  THE  MAXIMUM  RESIDUAL, 

SUM=O.ODOO 

DO  5  1*1, NSQ 

IF(0( I ).LE.O.ODOO)GO  TO  5 

SUM«SUM+Q( I ) 
5  CONTINUE 
C SET  PARAMETERS 

NSQ»N*N 

NSQP2N=NSQ+N+N 

NSQMN=NSO-N 

ITER=0 

ICNT=0 

WRITE(6,903) 
10  IF( ICNT.EO.M) ICNT=0 
--THE  ACTUAL  STONE  ITERATIVE  METHOD  BEGINS  HERE.  ICNT  KEEPS 
—TRACK  OF  THE  ITERATION  PARAMETERS.  FOR  A  DESCRIPTION  OF  THE 
--FUNCTIONS  OF  EACH  SUBROUTINE  AND  ITS  RELATION  TO  THE  ITERATIVE 
— METHODt  SEE  THE  COMMENTS  AT  THE  BEGINNING  OF  THE  INDIVIDUAL 
--SUBROUTINE. 

—EVEN  ITERATION  STEPS  ARE  PERFORMED  HERE. 

ICNT=ICNT+1 

ALPHA=PARAM( ICNT) 
—COMPUTE  THE  COEFFICIENTS  OF  THE  ALTERED  MATRIX,  A+N,  FROM  A. 

CALL  ALTERM(N,A1,A2,BL,CL,DL,EU,FU,ALPHA,NSQ,NSQP2N,NSQMN) 
—COMPUTE  THE  PRODUCT  A*T. 

CALL  PR0DMT(N,A1,A2,T,RN,NSQ,NSQP2N) 
—COMPUTE  THE  RESIDUAL  Q-A*T. 

CALL  RESIDL(RN,Q,NSQ) 
—COMPUTE  THE  MAXIMUM  NORMALIZED  RESIDUAL. 

CALL  MAXNR(NSQ,RN,SUM,BIG) 
—OUTPUT  THE  MAXIMUM  NORMALIZED  RESIDUAL  AFTER  ITER  ITERATIONS. 

WRITE<6,905> ITER, BIG 

ITER=ITER+1 
—SOLVE  L*V=RN  FOR  V. 

CALL  SOLVEL(N,BL,CL,DL,RN,NSQ,NSQMN) 
—SOLVE  U*DELTA=V  FOR  DELTA. 

CALL  SOLVEU(N,EU,FU,RN,NSQ,NSQMN) 
-COMPUTE  T(N+1)=T(N)+DELTA. 

CALL  NEWT(T,RN,NSQ> 
-TEST  FOR  END  OF  ITERATION  SCHEME 

IF( ITER.GE.ITMAX)  GO  TO  20 
-ODD  ITERATION  STEPS  ARE  PERFORMED  HERE.  THE  GRIDPOINTS  ARE 
—CONSIDERED  IN  A  NEW  ORDERING  FOR  THESE  STEPS(SEE  PAPER). 

-COMPUTE  THE  COEFFICIENTS  OF  THE  ALTERED  MATRIX,  A+N,  FROM  A. 

CALL  BALTRM(N,A1,A2,BL,CL,DL,EU,FU,ALPHA,NSQ,NSQP2N,NSQMN) 
—COMPUTE  THE  PRODUCT  A*T. 

CALL  PR0DMT(N,A1,A2,T,RN,NS0,NSQP2N) 
-COMPUTE  THE  RESIDUAL  0-A*T. 

CALL  RESIDL(RN,Q,NSQ) 
-COMPUTE  THE  MAXIMUM  NORMALIZED  RESIDUAL. 

CALL  MAXNR(NSO,RN,SUM,BIG) 
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C OUTPUT  THE  MAXIMUM  NORMALIZED  RESIDUAL  AFTER  ITER  ITERATIONS. 

WRITE(6,905)ITER,BIG 

ITER«ITER+1 
C SOLVE  L*V*RN  FOR  V. 

CALL  BSOLVL(N,BL»CLtDLtRN,NSQ,NSQMN) 
C SOLVE  U*DELTA=V  FOR  DELTA, 

CALL  BSOLVU(N,EUtFU,RNtNSOfNSOMN) 
C COMPUTE  T(N+1)»T(N)+DELTA. 

CALL  NEWT(T,RN,NSQ) 
C TEST  FOR  END  OF  ITERATION  SCHEME 

IF(  ITER.LT.ITMAX)GO  TO  10 

C COMPUTE  THE  PRODUCT  A*T. 

20  CALL  PRODMT(N»AltA2,T,RNrNSQ,NSQP2N) 
C COMPUTE  THE  RESIDUAL  Q-A*T. 

CALL  RESIDL(RN»Q,NSQ) 
C COMPUTE  THE  MAXIMUM  NORMALIZED  RESIDUAL  FOR  THE  LAST  ITERATION, 

CALL  MAXNR(NSOtRN,SUM,BIG) 

WRITE(6,905)ITER,BIG 
C OUTPUT  THE  FINAL  SOLUTION  VECTOR  T(I). 

WRITE(6,901) 

WRITE(6,902)  (  I  ,  T <  I  ) , I  =  1 ,NSQ ) 
C FORMAT  STATEMENTS 

901  FORMAT! 'OTHE  SOLUTION  VECTOR  T(I)  IS:') 

902  FORMAT (•  ',/,4('  T(',I3,»)  =  ', 016,8)) 

903  FORMAT( '0'tl3X, 'MAXIMUM  NORMALIZED  RESIDUAL') 
905  FORMAT('  I TERAT ION ' , I  3 1 6X, D16.8 ) 

STOP 
END 


SUBROUTINE  PARAMS ( N, Al , A2t M,P ARM,NSQP2N ) 

IMPLICIT  REAL*8( A-H,0-Z) 

DIMENSION  Al(NS0P2N)tA2(NS0P2N) »PARM(M) 
C 

C CALCULATE  THE  PARAMETERS  FOR  THE  USE  IN  ONE  CYCLE  OF  2*M 

C ITERATIONS  FOR  VARIABLE  Al  AND  A2  BY  TAKING  THE  ARITHMETIC 

C AVERAGE  OF  THE  MAXIMUM  PARAMETER  FOR  ALL  THE  GRID  POINTS, 

C 

PTMAX=0.0000 

NSQ=N*N 

DO  1  I=ltNSO 

IF( Al ( 1+1 >  ,EQ.0.0D00.OR.A2( I +N ) . EO .O.ODOO ) GO  TO  4 

TEMP1=A1 ( I+l )/A2( I+N) 

TEMP2=A2 ( I+N)/A1( I+l) 

GO  TO  3 
A  TEMP1=0.0D00 

TEMP2=0.0D00 
3  AMAX  =  1.0D00-DMIN1(  (2.0D00*( 1.0D00/(N+1 ) )  ) / ( 1 .0000+ ( TEMP2      )  )  « (2, 
X0D00*( 1.0000/(N+1  )  )  ) /( 1.0D00+( TEMPI      )  )  ) 

PTMAX=PTMAX+AMAX 
I  CONTINUE 

AMAX=PTMAX/DFLOAT(NSO) 

PARMI 1 )=0.0000 

H1=M-1 
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DO  2  1=2, M 

6X=DFLOAT< I-l)/DFLOAT(Ml) 
PARMU )*1.0DOO-( (1.0DOO-AMAX)**EX) 
2  CONTINUE 
M2=2*M 

WRITE(6,100)M2,(PARM( I ),I«1,M) 
100  FORMAT('0THE  UNORDERED  PARAMETERS  FOR  • 
XARE: • ,/,(D16.8)  ) 
RETURN 
ENO 


SUBROUTINE  ALT ERM( N ,A1 ,A2 ,BL ,CL,DL,EU, FU, ALPHA, NSQtNSQP2N ,NSQMN ) 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  AKNSQP2N)  ,  A2  (  NSQP2N  )  ,BL  (  NSOMN  )  ,CL  (  NSOMN  )  ,DL < NSO ) ,EU ( NSQ 
XMN) ,FU(NSOMN) 

-CALCULATE  THE  ALTERED  MATRIX  M+N=L*U  FROM  THE  GIVEN  MATRIX  M.  L  IS 
-A  TRI-DIAGONAL  MATRIX  CONSISTING  OF  THE  MAIN  DIAGONAL,  THE  LEFT 
-ADJACENT  DIAGONAL,  AND  A  LOWER  DIAGONAL  AT  A  DISTANCE  N  FROM  THE 
-MAIN  ONE.  THE  ELEMENTS  OF  THESE  DIAGONALS  ARE  STORED  IN  THE 
-VECTORS  DL,  CL,  AND  BL  RESPECTIVELY.  U  IS  A  TRI-DIAGONAL  MATRIX 
-CONSISTING  OF  THE  MAIM  DIAGONAL  EVERY  ELEMENT  OF  WHICH  IS  EQUAL 
-TO  THE  CONSTANT  ONE,  THE  RIGHT  ADJACENT  DIAGONAL,  AND  AN  UPPER 
-DIAGONAL  AT  A  DISTANCE  N  FROM  THE  MAIN  ONE.  THE  ELEMENTS  OF  THE 
-THE  TWO  OFF  DIAGONALS  ARE  STORED  IN  THE  VECTORS  EU  AND  FU 
-RESPECTIVELY.  VANISHING  ELEMENTS  OF  THE  DIAGONALS  ARE  NOT  STORED 
-IN  THESE  AFORE  MENTIONED  VECTORS. 

COMPUTATION  IS  DONE  IN  A  LEFT-TO-RIGHT  AND  DOWN-TO-UP  ORDER 
-SINCE  THE  ELEMENTS  OF  L*U  CALCULATED  AT  ANY  GRID  POINT  ONLY 
-DEPEND  ON  PREVIOUSLY  CALCULATED  ELEMENTS  AT  PREVIOUS  GRID  POINTS. 

PS,PW,PO,PE,  AND  PN  DEFINE  THE  LOWER  DIAGONAL,  LEFT  ADJACENT 
-DIAGONAL,  MAIN  DIAGONAL,  RIGHT  ADJACENT  DIAGONAL,  AND  UPPER 
-DIAGONAL  OF  THE  MATRIX  M. 

-DEFINE  CONSTANTS 

N1=N+1 

NN1=N1+N 

NN=N+N 

N21=NS0-N-N+1 

NSQM1=NSQ-1 
-GET  M+N=L*U  FROM  M 

PE  =  A1  (2I+AK3) 

PN=A2(N1 )+A2(NNl) 

PW=A1 (1)+A1<2) 

PS=A2(1)+A2(N1) 

PO=-(PE+PN+PW+PS) 
-THE  COMPUTATION  IS  PERFORMED  AT  THE  FIRST  POINT,  THEN  THE 
-INTERMEDIATE  POINTS,  AND  THEN  THE  LAST  POINT  OF  THE  GIVEN 
-HORIZONTAL  LINE  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  OF 
-COMPONENTS  OF  M  ALONG  THE  BOTTOM  AND  TOP  HORIZONTAL  LINES,  AND 
-THE  LEFT  AND  RIGHT  VERTICAL  LINES  OF  THE  GRID.  THE  ELEMENTS  OF  L*U 
-FOR  THE  FIRST  GRID  POINT  ARE  COMPUTED  HERE. 

0L( 1  )  =  P0 
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EU(1)=PE/DL(1) 

FU(1)=PN/DL(1) 

DO  15  1=2, N 

C 1  FOLLOWS  THE  REMAINING  GRID  POINTS  ON  THE  BOTTOM  HORIZONTAL 

C LINE.  THE  ELEMENTS  OF  L*U  FOR  THOSE  POINTS  ARE  COMPUTED  HERE.  I  IS 

C THE  BASIS  FOR  THE  INDICES  OF  Al,  A2,  AND  THE  ELEMENTS  OF  L*U. 

PE=AHI  +  l>+AlU  +  2) 

PN=A2< I*N)+A2< I+NN) 

PWsAl(I)+AHI  +  l) 

PS=A2(I )+A2U+N) 

PO=-(PE+PN+PW+PS) 

CLU-1 )«PW/Q.ODOO+ALPHA*FU(I-l>) 

DL(I  )=PO*CL(I-l)*(ALPHA*FU(I-l)-EU(I-l>) 

FUU )=(PN-ALPHA*CL( I-1)*FU< 1-11  )/DLU  ) 

C FOR  THE  LAST  GRID  POINT  ON  THE  BOTTOM  HORiZONTAL  LINE  THE  RIGHT 

C ADJACENT  DIAGONAL  ELEMENT  OF  L*U  VANISHES  SO  THE  COMPUTATION  OF 

C EU(N)  IS  INCORRECT  AND  WILL  BE  CORRECTED  IN  STATEMENT  16. 

EU(I )  =  PE/DL(  I) 

15  CONTINUE 

16  EU(N)=O.ODOO 

C L  AND  M  ARE  COUNTER  VARIABLES  THAT  AID  IN  CALCULATING  THE  CORRECT 

C SUBSCRIPTS  FOR  THE  COMPONENTS  OF  L*U.  M2  INCREMENTS  THE  Al  INDEX 

C AT  EACH  HORIZONTAL  LINE.  N1=N+1 ,N2 1=N**2-2*N+1 . 

M=0 

M2=0 

L  =  2 

DO  30  J=N1,N21,N 

C j  MOVES  THE  COMPUTATION  UP  TO  THE  NEXT  HORIZONTAL  LINE  OF  THE 

C GRID.  FOR  EACH  J  CALCULATIONS  ARE  DONE  AT  THE  FIRST  POINT,  THE 

C INTERMEDIATE  POINTS,  AND  THE  LAST  POINT  OF  THE  GIVEN  HORIZONTAL 

C LINE.  J  IS  THE  BASIS  FOR  THE  INDICES  OF  Al,  A2 ,  AND  THE  ELEMENTS 

C OF  L*U.  THE  CALCULATIONS  FOR  THE  FIRST  POINT  ON  EACH  HORIZONTAL 

C LINE  ARE  DONE  HERE. 

M2=M2+2 

PE=A1(J+M2+1)+A1( J+M2+2) 

PN=A2( J+N)+A2< J+NN) 

PW=A1(J+M2)+A1( J+M2+1) 

PS=A2( J)+A2< J+N) 

P0=-( PE+PN+PW+PS) 

BL(J-N)=PS/(1.0D00+ALPHA*EU(J-N-M) ) 

DL( J )=PO+BL( J-N)*( ALPHA*EU( J-N-M)-FU( J-N) ) 

EU(J-1-M)=(PE-ALPHA*BL( J-N)*EU( J-N-M) )/DL(J ) 

FU( J )=PN/DL( J) 

J1=J+1 

J2=J+N-2 

DO  25  I=J1,J2 

C 1  MOVES  THE  COMPUTATION  ALONG  THE  REMAINING  GRID  POINTS  ON  THE 

C GIVEN  HORIZONTAL  LINE  AND  THOSE  COMPUTATIONS  OF  THE  ELEMENTS  OF 

C L*U  APE  DONE  HERE.  I  IS  THE  BASIS  FOR  THE  INDICES  OF  Al,  A2,  AND 

C THE  ELEMENTS  OF  L*U. 

PE=A1(I+M2+1)+A1( I+M2+2) 

PN=A2( I+N)+A2( I+NN) 

PW  =  A1  (  I+M2  )+Al(  1+M2-H) 

K0=A2( I )+A2( I+N) 

PO=-(PE+PN+PW+PS) 

BL(I-N)-PS/( 1.0D00+ALPHA*EU( I-N-M) ) 

CL (  I-L  )=PW/<  1.0nOO+ALPHA*FU(  1-1 )  ) 
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DL(I  )*PO+BL(I-N>*<ALPHA*EU<I-N-M)-FU(I-N)  )+CL  (  I-L  )  *<  ALPHA*FU  (  1-1  )- 
XEUU-2-M) ) 

EU(I-1-M)*(PE-ALPHA*BL< 1 -N ) *EU ( I-N-M ) )/DL(I ) 

FU(  I  )  =  <PN-A|_PHA*CL(  I-L)*FU(I-1)  )/DL(I  ) 
25  CONTINUE 

I=J2+1 

THE  COMPUTATION  OF  THE  ELEMENTS  OF  L*U  FOR  THE  LAST  POINT  ON  THE 

GIVEN  HORIZONTAL  LINE  ARE  DONE  HERE  SEPARATELY  FROM  THE  OTHER  GRID 

POINTS  ON  THE  LINE  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING 

OF  THE  COMPONENT  EU  OF  L*U  AT  THIS  POINT. 

BL(J-1)=PS 

CL(  I-L)=PW/U.ODOO+ALPHA*FU(  1-1)  ) 

DLU  )  =  P0+BL<I-N)*(-FUU-N)>+CLU-L)*<ALPHA*FUU-i)-EU(I-2-M)  ) 

FU( I )=(PN-ALPHA*CL( I -L ) *FU ( I -1 ) ) /DL ( I ) 

M  =  M  +  1 

L  =  L  +  1  * 

30  CONTINUE 

1  =  1+1 

1    MOVES  THE  COMPUTATION  UP  TO  THE  FIRST  POINT  ON  THE  TOP 

HORIZONTAL  LINE  AND  THE  ELEMENTS  OF  L*U  FOR  THAT  POINT  ARE 

CALCULATED  HERE.  I  FORMS  THE  BASIS  FOR  THE  INDICES  OF  Alt  A2,  AND 

THE  ELEMENTS  OF  L*U. 

PE=A1( I+M2+1 )+AH I+M2+2) 

PN=A2( I+N)+A2( I+NN) 

PW=A1( I+M2)+A1( I+M2+1) 

PS  =  A2(I  )+A2(  I+N) 

PO=-(PE+PN+PW+PS) 

BL( I-N)=PS/( 1.0D00+ALPHA*EU( I-N-M) ) 

DL( I )=PO+BL( I-N)*( ALPHA*EU( I-N-M )-FU( I-N) ) 

EU( I-1-M)=(PE-ALPHA*BL( I-N ) *EU ( I -N-M ) )/DL( I ) 

11=1+1 

DO  45  K=I1,NSQM1 

K  MOVES  THE  COMPUTATION  ALONG  THE  REMAINING  GRID  POINTS  ON  THE 

TOP  HORIZONTAL  LINE  AND  THOSE  COMPUTATIONS  OF  THE  ELEMENTS  OF  L*U 

ARE  DONE  HERE.  K  FORMS  THE  BASIS  FOR  THE  INDICES  OF  Alt  A2t  AND 

THE  ELEMENTS  OF  L*U. 

PE  =  A1<K+M2  +  1  )+AKK  +  M2  +  2) 

PN=A2(K+N)+A2(K+NN) 

PW=A1 (K+M2 )+Al (K+M2+1) 

PS=A2(K)+A2(K+N) 

PO=-(PE+PN+PW+PS) 

8L(K-N)=PS/(1.0D00+ALPHA*EU(K-N-M) ) 

CL(K-N)=PW 

DL(K)=PO+BL(K-N)*(ALPHA*EU(K-N-M)-FU(K-N) ) +CL ( K-N ) * ( -EU ( K-N ) ) 

EU(K-N+1 )=(PE-ALPHA*BL(K-N)*EU(K-N-M) )/DL(K) 
45  CONTINUE 

THE  COMPUTATION  OF  THE  ELEMENTS  OF  L*U  FOR  THE  LAST  GRID  POINT 

OF  THE  TOP  HORIZONTAL  LINE  ARE  DONE  HERE  SEPARATELY  FROM  THE  OTHER 

GRID  POINTS  ON  THE  LINE  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE 

VANISHING  OF  THE  COMPONENTS  EU  AND  FU  OF  L*U  AT  THIS  POINT. 

BL(NSOMN)=PS 

CL(NSQMN)=PW 

DL(NSQ)=PO+BL(NSOMN)*(-FU(NSOMN) ) +CL ( NSQMN ) * (-EU ( NSQMN ) ) 

RETURN 

END 
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SUBROUTINE  PRODMT ( N, Al , A2, T ,RN,NSQ,NSQP2N ) 

IMPLICIT  REAL*8< A-H,0-Z) 

01  MENS  I  ON  A1(NSQP2N),A2(NSQP2N),T(NSQ),RN(NSQ> 
C 

C CALCULATE  THE  PRODUCT  M*T  WHERE  M  IS  A  MATRIX  AND  T  IS  A  VECTOR. 

C THE  COMPUTATION  IS  DONE  IN  A  LEFT-TO-RIGHT  AND  DOWN-TO-UP  ORDER 

C OVER  THE  GRID  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  OF 

C COMPONENTS  OF  M;  I.E.  THE  COMPUTATION  OF  RN»M*T  IS  PERFORMED 

C AT  THE  FIRST  GRID  POINT,  THEN  THE  INTERMEDIATE  GRID  POINTS,  AND 

C -THEN  THE  LAST  GRID  POINT  OF  EACH  HORIZONTAL  LINE  IN  ORDER  TO  TAKE 

C ADVANTAGE  OF  THE  VANISHING  OF  COMPONENTS  OF  M  ALONG  THE  BOTTOM 

C AND  TOP  HORIZONTAL  LINES  AND  THE  LEFT  AND  RIGHT  VERTICAL  LINES  OF 

C THE  GRID. 

C PS,PW,PO,PE,  AND  PN  DEFINE  THE  LOWER  DIAGONAL,  LEFT  ADJACENT 

C DIAGONAL,  MAIN  DIAGONAL,  RIGHT  ADJACENT  DIAGONAL,  AND  THE  UPPER 

C DIAGONAL  OF  THE  MATRIX  M. 

C 

C DEFINE  CONSTANTS 

N1=N+1 

NN1=N1+N 

NN=N+N 

N21=NSQ-N-N+1 
C BEGIN  PRODUCT  CALCULATION 

PE  =  A1  (2)+AK3) 

PN=A2(N1 )+A2(NNl) 

PW=A1(1)+A1(2) 

PS  =  A2U  )+A2(Nl) 

PO=-(PE+PN+PW+PS) 

C THE  CALCULATION  OF  RN=M*T  AT  THE  FIRST  GRID  POINT  IS  PERFORMED 

C HERE.  NOTE  THE  COMPONENT  PS  OF  M  VANISHES  ALONG  THE  BOTTOM 

C HORIZONTAL  LINE. 

RN(1 )=PO*T< 1)+PE*T(2)+PN*T(N1) 

DO  15  1=2, N 

c i    MOVES  THE  COMPUTATION  ALONG  THE  REMAINING  GRID  POINTS  ON  THE 

C BOTTOM  HORIZONTAL  LINE.  RN=M*T  AT  THESE  POINTS  IS  COMPUTED  HERE.  I 

c IS  THE  BASIS  FOR  THE  INDICES  OF  Al,  A2 ,  AND  RN. 

PE  =  A1  (  1+1  )+AH  1+2) 

PN=A2( I+N)+A2( I+NN) 

PW=A1 ( I )+Al( 1+1 ) 

PS=A2( I )+A2( I+N) 

PO=-(PE+PN+PW+PS) 
C FOR  I=N  THE  VALUE  OF  RN(N)  IS  INCORRECT  AND  IS  CORRECTED  BELOW. 

RN(I )=PW*T( 1-1 )+PO*T( I )+PE*T( I+1)+PN*T( I+N) 
15  CONTINUE 

C RN=M*T  AT  THE  LAST  GRID  POINT  ON  THE  BOTTOM  HORIZONTAL  LINE  IS 

C CALCULATED  HERE  SEPARATELY  FROM  THE  INTERMEDIATE  POINTS  ON  THE 

C LINE  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  COMPONENT  PE  OF 

C *.  THIS  IS  DOME  BY  SUBTRACTING  PE*T(N+1)  FROM  THE  RN(N)  JUST 

C CALCULATED. 

KN(N)=RN(N)-PE*T(N+1 ) 
C N1=N+1,  N21=N**2-2*N+1 

M1=N1 

M2  =  N 

30  J=N1,N21,N 
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c j  MOVES  THE  COMPUTATION  UP  TO  THE  NEXT  HORIZONTAL  LINE  OF  THE 

c GRID.  FOR  EACH  J  RN=M*T  IS  CALCULATED  AT  THE  FIRST  GRID  POINT,  THE 

c INTERMEDIATE  POINTS,  AND  THE  LAST  GRID  POINT  OF  THE  GIVEN 

c HORIZONTAL  LINE.  Ml,  M2  AND  J  ARE  THE  BASES  FOR  THE  INDICES  OF  Al, 

c A2,  AND  RN.  RN=M*T  FOR  THE  FIRST  GRID  POINT  OF  THE  HORIZONTAL  LINE 

c CALCULATED  HERE  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING 

C COMPONENT  PW  OF  M. 

MlxMl+3 

M2=M2+1 

PE=A1(M1)+A1(M1+1) 

PN=A2(M2)+A2(M2+N) 

PW=A1(M1-1)+A1(M1) 

PS=A2(M2-N)+A2(M2) 

PO=-(PE+PN+PW+PS) 

RN(J)=PS*T(  J-N)+PO*T(  J)+PE*T(J+l)«-PN*T(J  +  N> 

J1=J+1 

J2=J+N-1 

DO  25  I=J1,J2 

c 1  MOVES  THE  COMPUTATION  ALONG  THE  REMAINING  GRID  POINTS  ON  THE 

C GIVEN  HORIZONTAL  LINE  AND  RN=M*T  FOR  THOSE  POINTS  IS  COMPUTED 

C HERE.  Ml,  M2  AND  I  ARE  THE  BASES  FOR  THE  INDICES  OF  Al,  A2  AND  RN. 

M1=M1+1 

M2=M2+1 

PE  =  A1 (Ml )+Al (Ml  +  1  ) 

PN=A2(M2)+A2(M2+N) 

PW  =  A1  (M1-1)+A1(M1) 

PS=A2(M2-N)+A2(M2) 

PO=-(PE+PN+PW+PS) 

RN(I )=PS*T( I-N)+PW*T( I-l)+PO*T( I )+PE*T( 1+1 )+PN*T( I+N) 
25  CONTINUE 

C The  COMPUTATION  OF  RN=M*T  AT  THE  LAST  GRID  POINT  ON  THE  GIVEN 

c HORIZONTAL  LINE  IS  DONE  HERE  SEPARATELY  FROM  THE  OTHER  GRID  POINTS 

c ON  THE  LINE  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  COMPONENT 

C PE  OF  M.  THIS  IS  DONE  BY  SUBTRACTING  PE*T(J2+1)  FROM  THE  RN(J2) 

C JUST  CALCULATED. 

RN(J2)=RN( J2)-PE*T( J2+1) 
30  CONTINUE 

L=J2+1 

C l  MOVES  THE  COMPUTATION  UP  TO  THE  FIRST  GRID  POINT  ON  THE  TOP 

C HORIZONTAL  LINE  AND  RN  =  M*T  AT  THAT  POINT  IS  COMPUTED  HERE.  Ml,  M2 

C AND  L  ARE  THE  BASES  FOR  THE  INDICES  OF  Al,  A2  AND  RN.  NOTE  THE 

C COMPONENT  PN  OF  M  VANISHES  ALONG  THIS  LINE  AND  PW  VANISHES  AT 

C THIS  POINT. 

Ml=Ml+3 

M2=M2+1 

PE=A1 (Ml )+Al(Ml+l) 

PN=A2(M2)+A2(M2+N) 

PW=A1 (Ml-1 )+Al(Ml) 

PS=A2(M2-N)+A2(M2) 

P0=-( PE+PN+PS+PW ) 

RN(L )=PS*T(L-N)+P0*T(L)+PE*T(L+1 ) 

L1=L+1 

DO  45  K=L1,NS0 

C K  MOVES  THE  COMPUTATION  ALONG  THE  REMAINING  GRID  POINTS  ON  THE 

C TOP  HORIZONTAL  LINE  AND  RN=M*T  AT  THOSE  POINTS  IS  COMPUTED  HERE. 

C Ml,  M2  AND  K  ARE  THE  BASES  FOR  THE  INDICES  OF  Al,  A2  AND  RN. 

M1=M1+1 
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M2=M2+1 

PE=A1(M1)+A1(M1+1) 

PN=A2<M2)+A2(M2+N) 

PW=A1(M1-1)+A1(M1) 

PS=A2(M2-N)+A2(M2) 

PO=-(PE+PN+PW+PS) 
C FOR  K=NSQ  THE  VALUE  OF  RN(NSO)  IS  COMPUTED  BELOW. 

IF(K.EQ.NSQ)  GO  TO  40 

RN(K)=PS*T<K-N)+PW*T(K-1)+P0*T<K)+PE*T(K+1 ) 

GO  TO  45 

c THE  COMPUTATION  RN=M*T  AT  THE  LAST  GRID  POINT  ON  THE  TOP 

C HORIZONTAL  LINE  IS  PERFORMED  HERE  SEPARATELY  FROM  THE  OTHER  GRID 

c POINTS  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  OF  THE 

C COMPONENTS  PE  AND  PN  OF  M. 

40  RN(K)=PS*T(K-N)+PW*T(K-1)+P0*T(K) 
45  CONTINUE 

RETURN 

END 


SUBROUTINE  RES  I DL ( RN, Q, NSQ ) 

IMPLICIT  REAL*8( A-H,0-Z) 

DIMENSION  RN(NSO) tQ(NSQ) 
C 

c THIS  SUBROUTINE  CALCULATES  THE  RESIDUAL  RN=Q-M*T=Q-RN  WHERE  Q  IS 

c AN  INPUT  VECTOR  AND  RN=M*T  IS  A  VECTOR  COMPUTED  BY  THE  SUBROUTINE 

C PRODMT. 


C 


DO  1  I=ltNSO 
RN( I  )=Q( I )-RN( I ) 
1  CONTINUE 
RETURN 
END 


SUBROUTINE  MAXNR ( NSO t RN , SUM, B I G ) 

IMPLICIT  REAL*8( A-H,0-Z) 

DIMENSION  RN(NSQ) 
C 

C FIND  THE  MAXIMUM  NORMALIZED  RESIDUAL  AFTER  ITER  ITERATIONS  OF 

C THE  PROCEDURE. 

C 

BIG=O.ODOO 

DO    15    K  =  UNSO 

TEST=DABS(RN(K) ) 

IF(TEST.GT.BIG»BIG=TEST 
15    CONTINUE 

BIG-BIG/SUM 

RETURN 

FND 
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SUBROUTINE  SOLVEL  (  N,  BL  »CLt  DL,RNtNSQ  t  NSQMN  ) 

IMPLICIT  REAL*8( A-H,0-Z) 

DIMENSION  Bl < NS QMN ),CL< NSQMN ),DL(NSQ)tRN(NSQ) 

tutc  ciiRpnitTlNF  SOLVES  THE  EQUATION  L*V=RN  FOR  V  BY  DIRECTLY 
IIlIS^lSG^SE'^lTlS^F'ES^TlStis  GENERATED  BY  THE  ABOVE  EQUATION. 

L  IS  THE  TRI-DIAGONAL  MATRIX  CALCULATED  IN  SUBROUTINE  ALTERM  AND 

RN  IS  THE  VECTOR  CALCULATED  IN  SUBROUTINE  ^SIDL.  I.E. 

—  THE  FOLLOWING  RECURSION  FORMULA  IS  SOLVED  FOR  V  STARTING  AT  THE 

FIRST  GRID  POINT  AND  MOVING  IN  A  LEFT-TO-RIGHT  AND  DOWN-TO  UP 

ORDER  ACROSS  THE  GRID: 

BL(I )*V(I-N)+CL(I)*V(I-l)+OL(I )*V(I)=RN(I) 

—NOTE:  IN  ORDER  TO  SAVE  STORAGE  V  IS  RESTORED  IN  THE  VECTOR  RN. 
—ALSO,  REMEMBER  THAT  BL  AND  CL  ONLY  CONTAIN  NON-ZERO  VALUES  OF  THE 

—  DIAGONALS  SO  THEIR  SUBSCRIPTS  ARE  NOT  EXACTLY  THOSE  REPRESENTED 
IN  THE  ABOVE  RECURSION  FORMULA. 

RN(1)=RN(1)/DLU) 

DO  10  1=2, N 
RN(I)=(RN(I>-CL(I-1)*RN<I-1>>/DL(I) 

10  CONTINUE 

M  =  2 

L  =  N 

DO  30  1=2, N 
RN(L+1)=(RN(L+1)-BL(L+1-N)*RN(L+1-N>>/DL(L+1) 

RN(J°L)=(RN(J+L)-BL(J+L-N)*RN(J+L-N,-CL(J+L-M)*RN(J+L-1))/DL<J+L) 

20  CONTINUE 

M=M+1 

L  =  L+N 
30  CONTINUE 

RETURN 

END 


SUBROUTINE  SOL VEU ( N, EU, FU, V ,NSQ t NSQMN ) 

IMPLICIT  REAL*8( A-H,0-Z) 

DIMENSION  EU ( NSQMN ),FU( NSQMN) ,V(NSQ) 

C THIS  SUBROUTINE  SOLVES  THE  EQUATION  U*DELTA=V  FOR  DELTA  BY 

C_ DIRECTLY  SOLVING  THE  SYSTEM  OF  EQUATIONS  GENERATED  FROM  THE  ABOVE 

C EQUATION.  U  IS  THE  TRI-DIAGONAL  MATR I X  CALCULATED  I N  SUBROUT  NE 

C ALTERM  AND  V  IS  THE  VECTOR  CALCULATED  IN  SUBROUTINE  SOLVEL.  I.E. 

C THE  FOLLOWING  RECURSION  FORMULA  IS  SOLVED  FOR  DELTA  STARTING 

C AT  THE  LAST  GRID  POINT  AND  MOVING  IN  A  RIGHT-TO-LEFT  AND  UP  TO 

C DOWN  ORDER  ACROSS  THE  GRID: 

C DELTA ( I )+EU( I >*DELTA< 1+1 )+FU< I ) *DELTA ( I+N ) =V ( I ) 
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C 

C NOTE:  IN  ORDER  TO  SAVE  STORAGE  DELTA  IS  RESTORED  IN  THE  VECTOR  V 

C WHICH  IN  TURN  IS  RESTORED  IN  RN.  ALSOf  REMEMBER  THAT  EU  AND  FU 

C ONLY  CONTAIN  NON-ZERO  VALUES  OF  THE  DIAGONALS  SO  THEIR  SUBSCRIPTS 

C ARE  NOT  EXACTLY  THOSE  REPRESENTED  IN  THE  ABOVE  RECURSION  FORMULA. 

C 

NM1=N-1 

V(NSQ)=V(NSQ) 

DO  10  I=1»NM1 

N1=NSQ-I 

N3=NSQMN-I+1 

V(NI)=V(N1)-EU(N3)*V<N1+1) 
10  CONTINUE 

N1=NS0 

N2=NS0MN 

DO  30  1=1, NM1 

N1=N1-N 

V(N1 )=V(N1 )-FU(Nl)*V(Nl+N) 

DO  20  J=lfNMl 

N2=N2-1 

N3=N3-1 

V(N2)=V(N2 )-EU(N3)*V(N2+l)-FU(N2)*V(N2+N) 
20  CONTINUE 

N2=N2-1 
30  CONTINUE 

RETURN 

END 


SUBROUTINE  NEWT ( T  ,DELT A ,NSQ ) 

IMPLICIT  REAL*8(A-H,0-Z ) 

DIMENSION  T(NSQ> ,DELTA(NSQ) 
C 

C CALCULATE  THE  NEW  VALUE  OF  THE  VECTOR  T  FROM  THE  OLD  VALUE  AND 

C THE  NEW  INCREMENT,  DELTA .( REMEMBER  DELTA  IS  STORED  IN  V  WHICH 

C IN  TURN  IS  STORED  IN  RN)  I.E.  CALCULATE: 

C 

c T(N+1)=T(N)+DELTA(N+1)=T(N)+RN(N+1) 

C 

DO  1  1=1, NSO 

T( I )=T( I )+DELTA( I ) 
1  CONTINUE 

RETURN 

END 


SUBROUTINE  B ALTRM ( N , A  1 , A2 , B ,C , D, E , F , ALPHA , NS0,NSQP2N ,NSQMN ) 

IMPLICIT  RFAL*8(A-H,0-2 ) 

DIMENSION    AKNS0P2N) , A2 ( NS0P2N ) , B < NSQMN ) ,C(NSO) ,D(NSO) ,E(NSQ),F(NS 

/OMN) 
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C CALCULATE  THE  ALTERED  MATRIX  M+N»L*U  FROM  THE  GIVEN  MATRIX  M.  L  IS 

C A  TRI-DIAGONAL  MATRIX  CONSISTING  OF  THE  MAIN  DIAGONAL,  THE  LEFT 

C ADJACENT  DIAGONAL,  AND  A  LOWER  DIAGONAL  AT  A  DISTANCE  N  FROM  THE 

C MAIN  ONE.  THE  ELEMENTS  OF  THESE  DIAGONALS  ARE  STORED  IN  THE 

C VECTORS  D,  C,  AND  B  RESPECTFULLY.  U  IS  A  TRI-DIAGONAL  MATRIX 

C CONSISTING  OF  THE  MAIN  DIAGONAL  EVERY  ELEMENT  OF  WHICH  IS  EQUAL 

C TO  THE  CONSTANT  ONE,  THE  RIGHT  ADJACENT  DIAGONAL,  AND  AN  UPPER 

C DIAGONAL  AT  A  DISTANCE  N  FROM  THE  MAIN  ONE.  THE  ELEMENTS  OF  THE 

C the  TWO  OFF  DIAGONALS  ARE  STORED  IN  THE  VECTORS  E  AND  F 

C RESPECTFULLY.  VANISHING  ELEMENTS  OF  C  AND  E  ARE  STORED  AS  ZEROES 

C IN  THESE  AFORE  MENTIONED  VECTORS. 

c THE  COMPUTATION  IS  PERFORMED  IN  REVERSE  ORDER  FROM  SUBROUTINE 

C ALTERM.  THAT  IS,  THE  SOLUTION  IS  COMPUTED  AT  THE  GRID  POINTS 

c IN  THE  J=N  ROW  WITH  I  ASSUMING  THE  VALUES  1  UP  TO  N.  THEN 

C J  IS  DECREMENTED  BY  1  AND  I  AGAIN  ASSUMES  THE  RANGE  OF 

C VALUES  1  THROUGH  N.  THIS  LEFT-TO-RIGHT  AND  UP-TO-DOWN 

C ORDERING  OF  THE  GRID  POINTS  FOR  EVEN  ITERATIONS  BEGINNING  WITH 

C THE  SECOND  MAY  BE  CONSIDERED  A  REORDERING  OF  THE  ORIGINAL  POINTS. 

C PS,PW,PO,PE,  AND  PN  DEFINE  THE  LOWER  DIAGONAL,  LEFT  ADJACENT 

C DIAGONAL,  MAIN  DIAGONAL,  RIGHT  ADJACENT  DIAGONAL,  AND  UPPER 

C DIAGONAL  OF  THE  MATRIX  M. 

C 

C COMPUTE  CONSTANTS 

NM2=N-2 

NP1=N+1 

NSQP1=NSQ+1 

NSQPN=NSQ+N 

NM1S0=NSQ-N-N+1 

NSQMNl=NSOMN+l 
C THE  SUBSCRIPTS  OF  Al  AND  A2  FOLLOW  THE  NORMAL  ORDERING. 

PE=A1 (NSOPN)+Al (NSOPN+1) 

PN=A2(NSQP1)+A2(NSQP1+N) 

PW  =  A1  (NSQPN-1  )+AKNSQPN) 

PS=A2(NSOPl-N)+A2(NSOPl) 

PO=-(PE+PN+PW+PS) 

D(l )=PO 

E(1)=PE/D(1) 

F(l  )=PS/D(1) 

C(l  )=O.ODOO 

C OBSERVE  THAT  THE  MATRIX  M  IS  REORDERED  WITH  THE  RESULT  THAT  FOR 

C A  GIVEN  ROW  I,  PN  AND  PS  ARE  THE  COEFFICIENTS  OF  X(I-N)  AND  XU+N) 

C RESPECTFULLY.  THIS  REVERSES  THEIR  ROLE  IN  THE  OLD  ORDERING  WHERE 

C PN  AND  PS  ARE  THE  COEFFICIENTS  OF  XU+N)  AND  X(I-N). 

C THE  ROLE  OF  PE  AND  PW  IS  UNCHANGED. 

C THE  FOLLOWING  LOOP  COMPUTES  C,D,E,  AND  F  ALONG  THE  TOP 

C HORIZONTAL  LINE  OF  THE  GRID. 

DO  10  1=2, N 

C LOOP  PARAMETERS  FOLLOW  GRID  POINTS  IN  THE  NEW  ORDERING.  THESE 

C PARAMETERS  AVOID  SUPERFLUOUS  ADDITIONS  OF  I. 

IM1=I-1 

NS0I=NSQ+IM1 

NSQP1I=NSQP1+IM1 

NS0PNI=NSQPN+IM1 

PE=A1 (NSQPNI )+Al(NSQPNI+l) 

PN=A2(NSQP1I )+A2(NSQPlI+N) 

PW=A1 (NSOPNI-1 )+Al (NSQPNI ) 

PS=A2 (NSQP1I-N)+A2(NSQP1I ) 


56 

PO=-(PE+PN+PW+PS) 

C(I )=PW/(1.0D00+ALPHA*F( 1-1 ) ) 

D(I  )=P0+C(I)*(ALPHA*F(I-1>-E(I-1)> 

E(I )=PE/0(I ) 

C FOR  I=N  THIS  VALUE  OF  E  IS  INCORRECT  AND  IS  CORRECTED  IN 

C STATEMENT  11. 

FII  )  =  <PS-ALPHA*Cm*FU-l))/0(l> 

10  CONTINUE 

11  E(N)=O.ODOO 

C THE  NEXT  SEGMENT  COMPUTES  B,CtD,E,  AND  F  FROM  THE  SECOND 

C HORIZONTAL  LINE  (IN  UP-TO-DOWN  ORDER)  TO  AND  INCLUDING  THE  N-l 

C HORIZONTAL  LINE.  THE  SEGMENT  IS  A  NESTED  LOOP.  ONE  EXECUTION  OF 

C THE  OUTER  LOOP  COMPUTES  THESE  COMPONENTS  ALONG  THE  FIXED 

C HORIZONTAL  LINE.  CERTAIN  PARAMETERS  AND  COMPONENTS  ARE  ZERO  AT  THE 

C ENDPOINTS  OF  HORIZONTAL  LINES.  ACCORDINGLYt  THE  COMPUTATIONS 

C FOR  A  FIXED  HORIZONTAL  LINE  PROCEED  AS  FOLLOWS: 

c (D  c=o  AND  PW=o  AT  THE  LEFT  ENDPOINT.  THIS  FACT  ABOUT  PW 

C REQUIRES  THE  COMPUTATION  AT  THIS  POINT  TO  BE  A  SPECIAL  CASE 

C NOTE  THAT  SINCE  C=0  AT  THE  ENDPOINTS*  IT  IS  POSSIBLE  TO 

C AVOID  RESERVING  CORRESPONDING  STORAGE  IN  THE  C  ARRAY.  THIS 

C WAS  NOT  DONE  IN  ORDER  TO  SIMPLIFY  THE  CALCULATION  OF 

C SUBSCRIPTS  OF  C. 

C (2)  POINTS  2  THROUGH  N-l  ON  THE  GIVEN  HORIZONTAL  LINE  ARE 

C COMPUTED  ITERATIVELY. 

C (3)  E=0  AND  PE=0  AT  THE  RIGHT  ENDPOINTS.  PE=0  MAKES  THIS  A 

C SPECIAL  CASE.  NOTE  ALSO,  THAT  E=0  AT  THESE  ENDPOINTS  MAKES 

C IT  POSSIBLE  TO  AVOID  RESERVING  CORRESPONDING  STORAGE  IN  THE 

C E  ARRAY.  THIS  WAS  NOT  DONE  IN  ORDER  TO  SIMPLIFY  THE 

C CALCULATION  OF  THE  SUBSCRIPTS  OF  E. 

C (4)  NOTE  THAT  THE  VALUES  OF  C  AND  E  THAT  ARE  ZERO  AT  CERTAIN 

C ENDPOINTS  ARE  STORED  AS  ZERO  SO  THAT  IF  ONE  DESIRES  A 

C PRINTOUT  OF  THE  ARRAYS  FOR  DEBUGGING  PURPOSES  THE  CORRECT 

C VALUES  OF  C  AND  E  WILL  BE  STORED  FOR  THE  CORRESPONDING 

C ENDPOINTS. 

C THE  OUTER  LOOP  MOVES  THE  COMPUTATION  FROM  ONE  HORIZONTAL  LINE  TO 

C THE  NEXT  LOWER. 

Kl=-N-2 
C Kl  STEPS  Al  DOWN 

K2  =  -N 

C K2  STEPS  A2  DOWN.  THE  LOOP  VARIABLE  FOLLOWS  THE  INDEX  OF  POINTS 

C OF  THE  LEFT  VERTICAL  LINE  IN  THE  NEW  ORDERING. 

DO  30  I=NP1,NM1S0tN 

Kl=Kl+N+2 

K2=K2+N 

Ll=NS0-2-Kl 

C LI  IS  THE  INDEX  OF  Al  AT  THE  LEFT  GRID  POINT  OF  THE  CURRENT 

C HORIZONTAL  LINE. 

L2=NS0-N+1-K2 

C L2  IS  THE  INDEX  OF  A2  AT  THE  LEFT  GRID  POINT  OF  THE  CURRENT 

r HORIZONTAL  LINE. 

PE=A1 (LI )+Al (Ll+1 ) 

PN=A2(L2)+A2(L2+N) 

PW  =  A1  (Ll-1  )+AKLl  ) 

PS  =  A2(L2-N)+A2(L2  ) 

PO=-(PE+PN+PW+PS) 

[MN«I-N 

h( lMN)=PN/( 1 ,ODOO+ALPHA*E ( IMN) ) 
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C(I )=0.OD0O 

D(  I  )=PO+B( IMN)*( ALPHA*E( IMN)-F(IMN)) 

E( I )=<PE-ALPHA*B( I MN ) *E < I MN ) ) /D( I ) 

F(I )=PS/D(  I  ) 

c jj  FOLLOWS  THE  INDEX  OF  THE  INTERIOR  GRID  POINTS  IN  THE  NEW 

c ORDERING. 

DO  20  J=itNM2 

JJ=I+J 

L1J=L1+J 

L2J=L2+J 

PE  =  A1(L1J)+A1(LU  +  1) 

PN  =  A2(L2J)+A2U2J+N) 

PW=A1<L1J-1)+A1(L1J) 

PS=A2(L2J-N)+A2(L2J) 

PO=-(PE+PN+PW+PS) 

JMN=JJ-N 

B(JMN)=PN/(1.0D00+ALPHA*E< JMN) ) 

C( JJ)=PW/(1.0D00+ALPHA*F( JJ-1 ) ) 

D(JJ)=PO+B( JMN)*(ALPHA*E( JMN)-F(JMN) )+C ( J J ) * ( ALPHA*F ( J J-l )-E ( JJ-1 ) 
X) 

E(JJ  )=<PE-ALPHA*B( JMN)*E< JMN) )/D(JJ) 

F(JJ)=(PS-ALPHA*C( JJ)*F( JJ-1))/D(JJ) 
20  CONTINUE 

C THIS  SECTION  COMPUTES  THE  COMPONENTS  AT  THE  RIGHT  ENDPOINT  ON  THE 

C GIVEN  HORIZONTAL  LINE. 

L1J=L1J+1 

L2J=L2J+1 

PE=A1 (L1J)+A1(L1J+1) 

PN=A2(L2J )+A2(L2J+N) 

PW=A1(L1J-1)+A1(L1J) 

PS=A2(L2J-N)+A2(L2J) 

PO=-(PE+PN+PW+PS> 
C E  =  0  AND  PE  =  0  AT  THIS  ENDPOINT.  ZERO  IS  STORED  IN  E(JJ). 

JJ=JJ+1 

JMN=JJ-N 

B(  JMN)=PN 
C SINCE  E  =  0 

C(  J J  )=PW/(1.0D00  +  ALPHA*F( JJ-1 ) ) 

D(JJ)=PO+B(JMN)*(-F( JMN) ) +C ( J J ) * ( ALPHA*F ( J J-l )-E ( J J-l ) ) 

E(  JJ)=O.ODOO 
C SINCE  PW  DOES  NOT  APPEAR. 

F(JJ)=(PS-ALPHA*C< JJ)*F( JJ-1 ) )/D( JJ) 
30  CONTINUE 
C THE  COMPUTATION  IS  NOW  AT  THE  FINAL  HORIZONTAL  LINE. 

PE  =  A1  (2)+AK3) 

PN=A2(NP1 )+A2(NPl+N) 

PW  =  A1(1  )+AK2) 

PS=A2(NP1-N)+A2(NP1) 

PO=-(PE+PN+PW+PS) 
C NM1SQ=(N**2)-(2*N)  +  1 

B(NMlSQ)=PN/( 1.0D00+ALPHA*E(NM1S0) ) 
C NSQMN1  =  (N**2)-N+1 

D(NS0MN1 )=P0+B(NM1S0)*(ALPHA*E(NM1S0)-F(NM1S0) ) 

EINSQMN1 )=(PE-ALPHA*B(NM1S0)*E(NM1S0) )/D(NSOMNl) 

C(NS0MN1 )=0.0D00 

C THE  COMPUTATION  AT  THE  INTERIOR  POINTS  OF  THE  FINAL  HORIZONTAL 

C LINE  IS  D0NE  HERE.  THE  LOOP  VARIABLE  INDEXES  THE  GRID  POINTS 
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IN  THE  NEW  ORDERING.  K  ADJUSTS  THE  INDICES  Kl  AND  K2  OF  Al  AND  A2. 

K=-l 

L1=NSQMN1+1 

L2=NSQ-1 

DO  60  J=LltL2 

K  =  K  +  1 

Kl=3+K 

K2=N+2+K 

PE=A1(K1)+A1(K1+1> 

PN=A2(K2)+A2(K2+N) 

PW=A1 IK1-1)+A1(K1) 

PS=A2(K2-N)+A2(K2) 

PO=-(PE+PN+PW+PS) 

JMN=J-N  . 

B(JMN)=PN/(1.0D00+ALPHA*E( JMN)  ) 

C(J)=PW 

D(J)=PO+B( JMN)*( ALPHA*E(JMN)-F( JMN) ) +C ( J ) * ( -E ( J-l ) ) 

EU)  =  (PE-ALPHA*B<  JMN)*E(  JMN) )/D(J) 
60  CONTINUE 
WE  ARE  NOW  AT  THE  RIGHT  ENDPOINT.  NSOMN= ( N**2 )-N. 

K1=N+1 

K2=N+N 

PE  =  A1  (Kl  )+AKKl  +  l  ) 

PN=A2(K2)+A2(K2+N) 

PW=A1(K1-1)+A1(K1) 

PS=A2(K2-N)+A2(K2) 

PO=-(PE+PN+PW+PS> 

8(NSQMN)=PN 

C(NSO)=PW 

D(NSQ)=PO+B(NSQMN)*(-F(NSQMN) ) +C ( NSO ) * (-E ( NSO-1 ) ) 

E(NSQ)=O.ODOO 

RETURN 

END 


SUBROUTINE  BSOLVL ( N f B t C » D ,RN ,NSO t NSOMN ) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  B( NSOMN) ,C(NSQ) ,D(NSO) tRN(NSQ) 
C 

C THIS  SUBROUTINE  SOLVES  THE  EQUATION  L*V=RN  FOR  V  BY  DIRECTLY 

C SOLVING  THE  SYSTEM  OF  EQUATIONS  GENERATED  BY  THE  ABOVE  EQUATION, 

C 1_  is  THE  TRI-DIAGONAL  MATRIX  CALCULATED  IN  SUBROUTINE  BALTRM  AND 

C RN  IS  THE  VECTOR  CALCULATED  IN  SUBROUTINE  RESIDL.  I.E. 

C THE  FOLLOWING  RECURSION  FORMULA  IS  SOLVED  FOR  V  STARTING  AT  THE 

C FIRST  GRID  POINT  AND  MOVING  IN  A  LEFT-TO-RIGHT  AND  DOWN-TO-UP 

C ORDER  ACROSS  THE  GRID: 

C 

C BL(I)*V(I-N)+CL(I)*VU-l)+DLm*V(I)=RN(I) 

C 

C NOTE:  IN  ORDER  TO  SAVE  STORAGE  V  IS  RESTORED  IN  THE  VECTOR  RN. 

C 

fN(l)  =RN(  1  )/D(  1 ) 

DO  10  I=?,N 
Ml  I  )  =  (kN(  1  )-C(  I  )*RN(  1-1 )  )/D(  I  ) 
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10  CONTINUE 
N1=N+1 

00  20  I=N1,NSQ 
RN( I )=<RN( I )-B<I-N)*RN( I-N )-C ( I ) *RN ( 1-1 ) ) /D( I ) 

20  CONTINUE 
RETURN 
END 


SUBROUTINE  BSOL VU ( N, E t F , Vt NSO »NSQMN ) 

IMPLICIT  REAL*8< A-H,0-Z) 

DIMENSION  E(NSQ),F(NSQMN)tV(NSQ) 
C 

c this  SUBROUTINE  SOLVES  THE  EQUATION  U*DELTA=V  FOR  DELTA  BY 

C DIRECTLY  SOLVING  THE  SYSTEM  OF  EQUATIONS  GENERATED  FROM  THE  ABOVE 

c EQUATION.  U  IS  THE  TRI-DIAGONAL  MATRIX  CALCULATED  IN  SUBROUTINE 

C BALTRM  AND  V  IS  THE  VECTOR  CALCULATED  IN  SUBROUTINE  BSOLVL.  I.E. 

c the  FOLLOWING  RECURSION  FORMULA  IS  SOLVED  FOR  DELTA  STARTING 

C AT  THE  LAST  GRID  POINT  AND  MOVING  IN  A  RIGHT-TO-LEFT  AND  UP-TO- 

C DOWN  ORDER  ACROSS  THE  GRID: 

C 

C DELTA  (I  )+EU(  I)*DELTA(  I  +  D+FUII  )*DELTA<  I+N)=V(  I) 

C 

C NOTE:  IN  ORDER  TO  SAVE  STORAGE  DELTA  IS  RESTORED  IN  THE  VECTOR  V 

C WHICH  IN  TURN  IS  RESTORED  IN  RN. 

C 

NM1=N-1 

DO  10  I=1,NM1 

N1=NSQ-I 

V(N1 )=V(N1 )-E(Nl)*V(Nl+l ) 
10  CONTINUE 

NSQM1=NSQ-1 

DO  20  I=N,NSQM1 

N1=NSQ-I 

V(N1 )=V(N1)-E(N1)*V(N1+1)-F(N1>*V(N1+N) 
20  CONTINUE 

RETURN 

END 


APPENDIX  B 
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ADI  PROGRAM  (DIRICHLET; 
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C  ADI  PROGRAM 

c THIS  PROGRAM  COMPUTES  THE  SOLUTION  X  OF  THE  LINEAR  SYSTEM  OF 

c EQUATIONS  A*X  =  S  AS  THE  LIMIT  OF  THE  SEQUENCE  DEFINED  BY: 

c (H  +  W*I )*X(K+1/2)=S-(V-W*I )*X(K) 

c (V+W*I )*X(K+1)=S-(H-W*I  )*X(K+l/2) 

C 

IMPLICIT  REAL*8  <A-H,0-Z) 

DIMENSION  A1(960)»A2(960>,X(900),S<900),H<900),ABND<6) ,WH(33)» 

XWV(33) ,THETA(6) 
C INPUT  N,S»X» Alt A2tNIT,ABND»BBND,CBNDfDBND 

C 

c NOTE:  IT  IS  NECESSARY  FOR  THE  USER  TO  INPUT  THE  FOLLOWING 

C VARIABLES  INTO  THIS  PROGRAM: 

C N  =  NUMBER  OF  GRID  POINTS  ON  ONE  SIDE  OF  THE  SOLUTION  GRID 

C (THERE  ARE  N**2  GRID  POINTS  IN  THE  TOTAL  GRID  SYSTEM). 

C S  =  THE  VECTOR  CONTAINING  THE  RIGHT  HAND  SIDES  OF  THE  LINEAR 

C SYSTEM. 

C X  =  THE  INITIAL  VECTOR. 

c A1,A2  =  THE  VECTORS  OF  THE  DIFFERENTIAL  EQUATION  COEFFICIENTS 

c NIT  =  THE  POWER  OF  2  ITERATIONS  TO  BE  PERFORMED  BY  THIS 

C PROGRAM. 

C ABND(l)  =  UPPER  BOUND  ON  THE  EIGENVALUES  OF  H. 

C BBND  =  LOWER  BOUND  ON  THE  EIGENVALUES  OF  H. 

C CBND  =  UPPER  BOUND  ON  THE  EIGENVALUES  OF  V. 

C DBND  =  LOWER  BOUND  ON  THE  EIGENVALUES  OF  V. 

C 

c THE  INPUT  THAT  FOLLOWS  IS  A  SAMPLE  OF  THE  NECESSARY  DATA. 

N  =  30 

NIT  =  5 

NSQ=N*N 

NSQP2N=NSQ+N+N 

DN=DFLOAT(N) 

DNl=DFLOAT(N+l) 

PI=3. 1415926535 89793 

DO  200  I=1,NSQP2N 

Al ( I )=-.5D00 
200  A2( I )=-.5D00 

DO  201  I=lt900 

X( I )=0.0D00 
200  S( I )=1.0D00 

ABND(1 )=4.0D00*DSIN(DN*PI/(2.0D00*DN1 ) )**2 

BBND=4.0D00*DSIN(PI/(2.0D00*DN1 ) )**2 

CBND=ABND(1 ) 

DBND=BBND 

C COMPUTE  THE  SUM  OF  THE  POSITIVE  RIGHT  HAND  SIDES  FOR  USAGE  IN 

C THE  NORMALIZING  OF  THE  MAXIMUM  RESIDUAL. 

SUM=O.ODOO 

DO  3  1=1, NSQ 

IF(S(  I  )  .LE.O.ODOOGO  TO  3 

SUM=SUM+S( I ) 
3  CONTINUE 
C SET  PARAMETERS 

MCNT=0 

NSQ=N*N 

NSQMN=N*N-N 

NSQP2N=NSQ+N+N 
C THE  OPTIMAL  ITERATION  PARAMETERS  FOR  MIT  =  2**NIT  ITERATIONS  ARE 
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c COMPUTED  HERE  BY  SUBROUTINE  OPTW.  THEY  ARE  RETURNED  IN  THE 

C VECTORS  WH  AND  WV  RESPECTIVELY,  IN  DECREASING  ORDER  AND  ARE 

C APPLIED  TO  THE  ITERATION  STEPS  IN  THAT  SAME  ORDER, 

MIT=2**NIT 

MIT1=MIT+1 

NIT1=NIT+1 

CALL  0PTW(NIT1,MIT1,WV,WH,ABND,BBND,CBND,DBND,THETA,PMU) 

WRITE(6,102) 

c the  ACTUAL  ADI  ITERATION  BEGINS  HERE.  MCNT  KEEPS  TRACK  OF  THE 

C ITERATION  PARAMETERS.  FOR  A  DESCRIPTION  OF  THE  FUNCTIONS  OF  EACH 

C SUBROUTINE  AND  ITS  RELATION  TO  THE  ITERATIVE  PROCEDURE,  SEE 

C the  COMMENTS  AT  THE  BEGINNING  OF  THE  INDIVIDUAL  SUBROUTINES. 

4  MCNT=MCNT+1 
c THE  FIRST  HALF  STEP  OF  THE  ITERATIVE  SEQUENCE  IS  PERFORMED  HERE. 

W=WH(MCNT) 
C COMPUTE  THE  PRODUCT  S-( V-W*I ) *X ( K ) . 

CALL  PROD V(N,S, NSQ, A2,W,X,NSQP2N) 
C COMPUTE  THE  INVERSE  OF  (H+W*I)  BY  TR I  ANGULAR  I ZAT ION . 

CALL  HH1 (N,A1,NSQMN,W,H,NSQ,NSQP2N) 
C COMPUTE  X(K+l/2)  BY  BACKWARD  SUBSTITUTION. 

CALL  PPKN,A1,NSQMN,W,H,NSQ,X,NSQP2N) 

C 

c THE  SECOND  HALF  STEP  OF  THE  ITERATIVE  SEQUENCE  IS  PERFORMED  HERE. 

W=WV(MCNT) 
C COMPUTE  THE  PRODUCT  S- (H-W*I ) *X ( K+l /2 ) . 

CALL  PRODH(N,S,NSQ,Al,NSQMN,W,X,NSQP2N) 
c COMPUTE  THE  INVERSE  OF  <V+W*I)  BY  TR I  ANGULAR  I ZAT ION. 

CALL  HH2(N,A2,NSQMN,W,H,NSQ,NSQP2N) 
c COMPUTE  X(K  +  1)  BY  BACKWARD  SUBSTITUTION. 

CALL  PP2(N,A2,NSQMN,W,H,NSQ,X,NSQP2N) 

c COMPUTE  THE  MAXIMUM  NORMALIZED  RESIDUAL  AND  THE  LI  VECTOR 

c NORM  AFTER  EACH  FULL  ITERATION. 

CALL  PRODMT(N,Al,A2,X,H,NSQ,NSQP2N) 

CALL  RESIDL!H,S,NSQ) 

BIG=O.ODOO 

XNORM=O.ODOO 

DO  5  1=1, NSQ 

TEMP  =  DABS(H(  I  )  ) 

XNORM=XNORM+TEMP 

IF(TEMP.GT.BIG)BIG=TEMP 
5  CONTINUE 

BIG=BIG/SUM 

WRITE(6,103)MCNT,BIG,XN0RM 
C TEST  CONVERGENCE  CRITERIA 

IF(MCNT.LT.MIT)  GO  TO  4 
C OUTPUT  THE  FINAL  SOLUTION  VECTOR  X(I). 

WRITE(6,100) 

WRITE(6,101)  ( I,X( I ) ,1=1, NSQ) 

C FORMAT  STATEMENTS 

LOO  FORMAT('OTHE  SOLUTION  VECTOR  X(I)  IS:'/) 

101  FORMAT!'  ',/,<►!'  X(M3t'l  =  ',016.8)) 

102  FORMAT!  '0',13X, 'MAXIMUM  NORMALIZED  RESI DUAL •  ,  3X  ,  • L  1  VECTOR  NORM') 

103  FORMAT!'  ITERATION' , 13, 6X,D16. 8) 
STOP 

^NO 


63 

SUBROUTINE  OPTW ( NI Tl »MI Tl , WVt WH t A, BTC* D» THETA, PMU ) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  WV(MIT1)*WH(M!T1) » A ( NI Tl ) , THETA ( NI Tl ) 

c this  SUBROUTINE  COMPUTES  THE  OPTIMAL  PARAMETERS  OF  THE  H  AND  V 

C MATRICES  FOR  USAGE  BY  THE  A.D.I.  PROGRAM,  THE  COMPUTATION  IS 

C BASED  ON  THE  NUMBER  OF  ITERATIONS*  MIT=2**NIT,  AND  IS  PERFORMED 

C ACCORDING  TO  THE  METHOD  DERIVED  BY  WACHSPRESS. 

NIT=NIT1-1 

MIT=MIT1-1 
C COMPUTE  BRACKET  RECURSION  RELATIONSHIPS 

THETA (1 )=(A(1 )*B-C*D)/(A( D+B+C  +  D) 

AP=A(1 )-THETA(l) 

BP=B-THETA(1) 

CP=C+THETA(1) 

DP=D+THETA(1) 

DO  1  I=2*NIT1 

A( I  )=DSQRT(AP*BP) 

B=(AP+BP)/2.0D00 

C=DSORT (CP*DP) 

D=(CP+DP)/2.0DOO 

THETA ( I  )  =  (A(  I  )*(B-D) )/(2.0*A( I )+B+D) 

AP=A( I )-THETA( I) 

BP=B-THETA( I ) 

CP=C+THETA( I ) 

DP=D+THETA( I ) 
1  CONTINUE 
C COMPUTE  ITERATION  PARAMETER  RECURSION  RELATIONSHIPS 

WP=DSQRT(AP*BP) 

WH(1 )=WP-THETA(NIT1 ) 

WV(1 )=WP+THETA(NIT1) 

NN=1 

DO  6  I=2»NIT1 

NN=NN*2 

NNN=NN+1 

K=NITl-I+2 

KK=K-1 

DO  3  INC=2»NNt2 

J=NN-INC+2 

JJ=J/2 

TEMP=WH( JJ)**2-A(K)**2 

IF(TEMP.LT,O.ODOO)TEMP=DABS(TEMP) 

WH(J)=WH( JJ)+DSQRT(TEMP) 

TEMP=WV( J J )**2-A(K)**2 

IF(TEMP.LT.O.ODOO)TEMP=DABS<TEMP) 

WV(J )=WV(JJ)+DSORT(TEMP) 

3  CONTINUE 

DO  4  J=3*NNNt2 
WH(J)=(A(K)**2)/WH(J-1) 
WV(J)=(A(K)**2)/WV( J-ll 

4  CONTINUE 

DO  5  J=1,NN 

WH( J)=WH( J+l) -THETA (KK) 

WV(J )=WV(J+1 )+THETA(KK) 

5  CONTINUE 
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6  CONTINUE 
ARRANGE  THE  OPTIMAL  PARAMETERS  IN  DECREASING  ORDER 

CALL  DECR(NN,WH,MIT1) 

WRITE(6»101)  <WH( I),I=1,NN) 

CALL  DECR(NN,WV,MIT1) 

WRITE(6,102)  (WVU)fl-ltNN)   . 
COMPUTE  THE  SPECTRAL  RADIUS 

PMU=( (WP-AP)/(WP+AP)>*< <WP-CP)/(WP*CP>> 

101  FORMATCOTHE  OPTIMAL  PARAMETERS  WH  IN  DECREASING  ORDER  ARE:'/,(D16 

X    8  )  ) 

102  FORMATCOTHE    OPTIMAL    PARAMETERS    WV    IN    DECREASING    ORDER    AREC/,(D16 

X.8)  ) 

103  FORMATCOTHE  SPECTRAL  RADIUS  =«,D16.8) 
RETURN 

END 


SUBROUTINE  DECR ( NN  ,WH , MIT1 > 

IMPLICIT  REAL*8  <A-H,0-Z) 

DIMENSION  WH(MITl) 

c THIS  SUBROUTINE  PLACES  THE  OPTIMAL  PARAMETERS  FOR  THE  H  AND  V 

C MATRICES  INTO  DECREASING  NUMERICAL  ORDER  FOR  USAGE  BY  THE  A.D.I 

C PROGRAM. 

NNM1=NN-1 

DO  1  I=ltNNMl 

11=1+1 

00  1  J=I1»NN 

IF(WHU  ).GE.WH(  J)  )G0  TO  1 

TEMP=WH(I ) 

WH( I  )=WH( J) 

WH(J)=TEMP 
1  CONTINUE 

RETURN 

END 


SUBROUTINE  PRODV < N ♦ S ,NSQ ♦ A2 »W , X , NS0P2N > 
IMPLICIT  REAL*8(A-H,0-Z) 
DIMENSION  S(NSO) , X ( NSO ) t A2 ( NS0P2N ) 

C 

c CALCULATE  THE  PRODUCT  ( S- < V-W*I ) *X ) .  COMPUTATION  OF  THE 

c COMPONENTS  IS  DONE  IN  DOWN-TO-UP  AND  LEFT-TO-RIGHT  ORDER  TO 

C TAKE  ADVANTAGE  OF  THE  VANISHING  OF  COMPONENTS  OF  V  ALONG  THE 

C BOTTOM  OR  TOP  LINES  OF  THE  GRID. 

C PS,PN,AND  PO  DEFINE  THE  LOWER  DIAGONAL,  UPPER  DIAGONAL*  AND 

c HAiN  DIAGONAL  ( SOUTHfNORTHt ORIGIN)  OF  V. 

C 

DO  2  1=1. N 

r i  FOLLOWS  THE  LOWER  HORIZONTAL  LINE.  FOR  EACH  I,  THE  COMPONENTS 

rn^  PRODUCT  ARE  COMPUTED  ALONG  THE  CORRESPONDING  VERTICAL  LINE. 
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c THE  FIRST  COMPUTATION  FOR  EACH  It  ANO  THE  LAST,  ARE  SPECIAL 

C CASES  OCCURRING  WHEN  PS  ANO  PN  ARE  ZERO. 

c I  MOVES  THE  COMPUTATIONS  FROM  LEFT  TO  RIGHT  ALONG  THE 

C (BOTTOM)  HORIZONTAL  LINE. 

IPN=I+N 

PS=A2(I  )+A2(IPN) 

PN=A2(IPN)+A2( IPN+N) 

PO=-(PS+PN) 
C SAVE  X(I)  TO  USE  AGAIN 

X1  =  X(I  ) 

X(I )«S(I )-< (PO-W)*X(I)+PN*X( IPN)) 

c THE  VARIABLE  OF  THIS  LOOP  STEPS  UP  EACH  VERTICAL  LINE  IN 

C ORDER  FROM  LEFT  TO  RIGHT.  THE  LOOP  VARIABLE  IS  THE  INOEX  OF  X. 

C HENCE  IPN  IS  THE  INDEX  OF  X  AT  THE  LOWEST  POINT  AND  LI  IS  THE 

C INDEX  OF  X  AT  THE  HIGHEST  POINT. 

L1=NSQ-N-N+I 

DO  1  J=IPN,LltN 

PS  =  PN 
C JPN  IS  THE  INDEX  OF  A2 

JPN=J+N 

PN=A2(JPN)+A2( JPN+N) 

PO=-(PS+PN) 

X2=X(J) 

X(J)=S(  J)-(PS*Xl+(PO-W)*X(  J)+PN*X<  JPN)  ) 

X1  =  X2 

1  CONTINUE 

C THIS  SEGMENT  COMPUTES  THE  COMPONENTS  ALONG  THE  UPPER 

C HORIZONTAL  LINE. 

PS  =  PN 

C K2  IS  THE  X  INDEX 

C |_1  =  N**2-2*N+K 

K2=L1+N 
C K3  IS  THE  A  INDEX 

K3=K2+N 

PN=A2(K3)+A2(K3+N) 

PO=-(PS+PN) 

X(K2)=S(K2)-(PS*Xl+(PO-W)*X(K2) ) 

2  CONTINUE 
RETURN 
END 


SUBROUTINE  HH1 ( N, Al , NSQMN, W ,H ,NSQ , NSQP2N ) 

IMPLICIT  REAL*8( A-HtO-Z) 

DIMENSION  H(NSQ)»A1(NSQP2N) 
C 

C THIS  SUBROUTINE  COMPUTES  THE  SOLUTION  OF  <H+W*I)X=S  BY  GAUSSIAN 

C ELIMINATION.  THE  MATRIX  <H+W*I)  IS  TRIDIAGONAL  WITH  THE  MAIN 

C DIAGONAL  AND  THE  TWO  ADJACENT  DIAGONALS  ON  EITHER  SIDE  OF 

C THE  MAIN  DIAGONAL. 

C COMPUTATION  IS  DONE  IN  THE  FOLLOWING  MANNER:  FIRST  FORWARD 

C ELIMINATION  IS  PERFORMED  TRANSFORMING  THE  TRIDIAGONAL  MATRIX 

C INTO  AN  UPPER  DIAGONAL  MATRIX  (SUBROUTINE  HH1 ) .  THEN  (SUBROUTINE 

C ppD  THE  RIGHT  HAND  SIDE  VECTOR  IS  TRANSFORMED,  THEN  THE 


66 

C BACKWARD  SUBSTITUTION  IS  PERFORMED  TO  SOLVE  THE  RESULTING 

c UPPER  TRIANGULAR  SYSTEM.  H  IS  SET  EQUAL  TO  ZERO  ALONG  THE 

C RIGHTMOST  VERTICAL  LINE  OF  THE  GRID,  WHICH  RESULTS  FROM  THE 

C FACT  THAT  THE  RIGHT  ADJACENT  DIAGONAL  VANISHES  AT  THESE  POINTS. 

C 

DO  1  K»N,NSQ,N 

H(K)=0.0 

1  CONTINUE 

c NExT  THE  COMPUTATION  IS  PERFORMED  ON  THE  FIRST  VERTICAL  LINE  OF 

C THE  GRID  SINCE  H  AT  THESE  POINTS  DEPENDS  ONLY  ON  THE  DIAGONAL 

C ELEMENT  AND  THE  RIGHT  ADJACENT  DIAGONAL  ELEMENT,  AND  DOES  NOT 

C DEPEND  QN  ANY  PREVIOUS  H. 

c PW,PE,AND  PO  DEFINE  THE  LEFT  DIAGONAL,  RIGHT  DIAGONAL,  AND 

C MAIN  DIAGONAL  ELEMENTS  OF  THE  MATRIX  H. 

NTOP=NSQMN+l 

M  =  l 

DO  2  I=l,NTOP,N 

C I  IS  THE  H  INDEX  AND  MOVES  THE  COMPUTATION  UP  THE  FIRST  VERTICAL 

C LINE.  IM  IS  THE  Al  INDEX. 

IM=I+M 

PW=A1( IM-1 )+Al ( IM) 

PE=A1(IM)+A1( IM+1) 

PO=-(PE+PW) 

H( I )=-PE/(PO+W) 

M=M  +  2 

2  CONTINUE 

C THE  COMPUTATION  OF  H  IS  NOW  PERFORMED  IN  LEFT-TO-RIGHT  AND  DOWN- 

C TO-UP  ORDER;  I.E.  H  IS  CALCULATED  ALONG  THE  REMAINING  GRID  POINTS 

C ON  THE  FIRST  HORIZONTAL  LINE  STARTING  WITH  THE  FIRST  AND  I  MOVES 

C THE  COMPUTATION  UP  TO  THE  NEXT  HORIZONTAL  LINE. 

Nl=N+l 

M  =  l 

11=2 

INM1=N-1 

DO  5  1=1, N 

PW=A1(I+M)+A1<  I+M+l) 

DO  4  J=II, INM1 
C J  IS  THE  H  INDEX  AND  Jl  IS  THE  Al  INDEX 

J1=J+1 

PE  =  A1  UD+AlUl  +  l) 

PO=-(PE+PW) 

H( J )=-PE/(PO+W+PW*H( J-l) > 

PW  =  PE 
4-  CONTINUE 

I  1  =  1 I+N 

INM1=INM1+N 

M=M+N1 
5  CONTINUE 

RETURN 

END 


MJBROUT INE    PP1 ( N , A  1 , N SOMN , W ,H ,NSO , X ,NSQP2N ) 
IMPLICIT    REAL*8 (A-HtO-Z» 
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DIMENSION  H(NSO) , X ( NSO > , Al ( NSOP2N ) 

C 

c this  SUBROUTINE  CONTINUES  THE  GAUSSIAN  ELIMINATION  PROCESS 

c BY  CALCULATING  THE  RECURSION  FORMULA  FOR  THE  RIGHT  HAND  SIDES; 

c i.e.  it  CALCULATES  P(M)=X(M)  FOR  M«1,...,N  BY  THE  FOLLOWING 

C FORMULA: 

C 

C P(M)*X(M)«(X(M)-PW*X(M-l))/(PO*W+PW*H(M-n ) 

C WHERE 

C P(l)»X(l)=X(l)/(PO^W) 

c 

C AND  WHERE  PW,  PE,  AND  PO  ARE  THE  LEFT  DIAGONAL,  RIGHT  DIAGONAL, 

C AND  MAIN  DIAGONAL  ELEMENTS  OF  H  AT  POINT  M. 

C COMPUTATION  IS  DONE  IN  A  LEFT-TO-RIGHT  AND  DOWN-TO-UP  ORDER 

C 0vER  THE  GRID.  FOR  EACH  HORIZONTAL  LINE  P(M)=X(M)  IS  CALCULATED 

C AT  THE  FIRST  POINT  AND  THEN  AT  THE  REMAINING  POINTS  IN  ORDER  TO 

c TAKE  ADVANTAGE  OF  THE  FACT  THAT  P(M)  DOES  NOT  DEPEND  ON  A  PREVIOUS 

C H  AT  THE  FIRST  POINT  AND  DOES  DEPEND  ON  A  PREVIOUS  H  AT  THE 

C REMAINING  POINTS.  THE  I  LOOP  MOVES  THE  COMPUTATION  UP  TO 

C THE  NEXT  HORIZONTAL  LINE  OF  THE  GRID, 

C 

INN=1 

11=2 

IN  =  N 

M  =  l 

DO  2  1=1, N 

C COMPUTATION  FOR  THE  FIRST  POINT  ON  THE  GIVEN  HORIZONTAL  LINE 

C is  PERFORMED  HERE.  INN  IS  THE  X  INDEX  AND  IM  IS  THE  Al  INDEX. 

IM=INN+M 

PW  =  A1(  IM-1  )+AK  IM) 

PE=A1(IM)+A1( IM+1) 

PO=-(PE+PW) 

X( INN)=X( INN)/(PO+W) 

DO  1  K=II,IN 

C THIS  LOOP  PERFORMS  THE  COMPUTATION  AT  THE  REMAINING  GRID  POINTS 

C ALONG  THE  GIVEN  HORIZONTAL  LINE.  K  IS  THE  X  INDEX  AND  KM  IS  THE 

C Al  INDEX. 

KM=K+M 

PW  =  PE 

PE=A1(KM)+A1(KM+1) 

PO=-(PE+PW) 

KK=K-1 

X(K)=(X(K)-PW*X(KK)  )/{PO+W  +  PW*H(KK)  ) 

1  CONTINUE 
INN=INN+N 
II=II+N 
IN=IN+N 
M=M+2 

2  CONTINUE 

C THIS  SECTION  COMPUTES  THE  BACKWARD  SUBSTITUTION  FOR  THE  GAUSSIAN 

C ELIMINATION  PROCEDURE;  I.E.  IT  COMPUTES: 

C 

C X(N)=P(N)=X(N) 

C X(M)=P<M)+H(M)*X(M+1)=X(M)+H(M)*X(M+1 )  FOR  M=N-1,...,1 

C 

C XI  SAVES  X(NSO)  TO  USE  AGAIN 

X1=X(NS0) 
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I1=NSQ-1 
00  3  1*1*11 

•KK  IS  THE  DECREASING  INDEX  FOR  X  AND  H 
KK=NSO-I 

X(KK)=X(KK)+H(KK)*X1 
-XI  SAVES  X(KK)  TO  USE  AGAIN 
X1«X(KK) 
CONTINUE 
RETURN 
END 


SUBROUTINE  PRODH ( N ,S ,NSO t Al ,NSQMN»W,X,NSQP2N ) 
IMPLICIT  REAL*8(A-H,0-Z) 
DIMENSION  S(NSQ),X(NSQ),A1(NSQP2N) 
C 

C CALCULATE  THE  PRODUCT  ( S-(H-W*I ) *X ) .  COMPUTATION  OF  THE  COMPONENTS 

C IS  DONE  IN  LEFT-TO-RIGHT  AND  DOWN-TO-UP  ORDER  TO  TAKE  ADVANTAGE 

C OF  THE  VANISHING  OF  COMPONENTS  OF  H  ALONG  THE  LEFTMOST  AND 

C RIGHTMOST  VERTICAL  LINES  OF  THE  GRID. 

C PW,PE,AND  PO  DEFINE  THE  LEFT  DIAGONAL*  RIGHT  DIAGONAL*  AND  MAIN 

C DIAGONAL  ( WEST, EAST , OR IGEN )  OF  H. 

C 

M=l 

NT0P=NSQMN+1 

DO  2  I=ltNTOPtN 

C i  FOLLOWS  THE  LEFTMOST  VERTICAL  LINE  OF  THE  GRID.  FOR  EACH  I  THE 

C COMPONENTS  OF  THE  PRODUCT  ARE  COMPUTED  ALONG  THE  CORRESPONDING 

C HORIZONTAL  LINE. 

C THE  FIRST  COMPUTATION  FOR  EACH  I*  AND  THE  LAST  ARE  SPECIAL 

C CASES  OCCURRING  WHEN  PW  AND  PE  ARE  ZERO.  THIS  FIRST  SECTION 

C PERFORMS  THE  COMPUTATION  AT  THE  FIRST  POINT  ON  THE  GIVEN 

C HORIZONTAL  LINE. 

C I  IS  jhe  X  INDEX  AND  IM  IS  THE  Al  INDEX.  XI  SAVES  X(I)  TO  USE 

C AGAIN 

IM=I+M 

X1=X( I ) 

PW=A1( IM-1 )+Al ( IM) 

PE=A1 ( IM)+A1< IM+1) 

PO=-(PE+PW) 

X(  I  )=S( I  )-(  (PO-W)*X( I )+PE*X( 1  +  1) ) 

C THIS  SECTION  PERFORMS  THE  COMPUTATION  ALONG  THE  INTERMEDIATE 

C POINTS  ON  THE  GIVEN  HORIZONTAL  LINE. 

11=1+1 

I2=I+N-2 

DO  1  J=Ilt 12 

C J  IS  THE  X  INDEX  AND  JM  IS  THE  Al  INDEX.  X2  SAVES  X(J)  TO  USE 

C AGAIN. 

JM=J+M 

X2  =  X(J  ) 

PW  =  PF 

PE=A1 ( JM)+A1 ( JM+1 ) 

PO=-(PE+PW) 

/  (J  )*S< J)-(PW*X1+(P0-W)*X< J)+PE*X(J  +  1 ) ) 
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X1  =  X2 

1  CONTINUE 

C THIS  SECTION  PERFORMS  THE  COMPUTATION  AT  THE  LAST  POINT  OF  THE 

C GIVEN  HORIZONTAL  LINE. 

L=I2+1 

LM=L*M 
C L  IS  THE  X  INDEX  AND  LM  IS  THE  Al  INDEX. 

PW=PE 

PE=A1(LM)+A1(LM*1> 

PO=-(PE+PW) 

X(L)=S(L)-(PW*XH-(PO-W)*X(L)  ) 

M  =  M  +  2 

2  CONTINUE 
RETURN 
END 


SUBROUTINE  HH2 ( N , A2 ,NSQMN , W ,H ,NSQ,NSQP2N ) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  H(NSQMN) ,A2(NSQP2N) 
C 

C THIS  SUBROUTINE  SOLVES  (V+W*I)X=S  BY  GAUSSIAN 

C ELIMINATION.  SINCE  V  IS  A  TRI-DIAGONAL  MATRIX  CONSISTING  OF  A 

C MAIN  DIAGONAL  AND  TWO  OFF  DIAGONALS  AT  A  DISTANCE  N  FROM  THE 

C MAIN  DIAGONAL,  GAUSSIAN  ELIMINATION  CAN  BE  PERFORMED  DIRECTLY 

C ON  THE  MATRIX. 

C FIRST  THE  FORWARD  ELIMINATION  PROCESS  IS  PERFORMED 

C TRANSFORMING  THE  TRI-DIAGONAL  MATRIX  INTO  AN  UPPER  DIAGONAL  MATRIX 

C (SUBROUTINE  HH2 )  THEN  (SUBROUTINE  PP2)  THE  RIGHT  HAND  SIDE  VECTOR 

C IS  TRANSFORMED  IN  ACCORDANCE  WITH  THE  UPPER  DIAGONAL  MATRIX;  I.E. 

C THE  P(I)=X(I)  ARE  CALCULATED.  THEN  THE  BACKWARD  SUBSTITUTION 

C (SUBROUTINE  PP1)  IS  PERFORMED  TO  GET  THE  RESULT. 

C ALONG  THE  FIRST  HORIZONTAL  LINE  OF  THE  GRID  THE  VALUES  OF 

C THE  UPPER  DIAGONAL  COMPONENT,  H,  DO  NOT  DEPEND  ON  ANY  PREVIOUS 

C VALUES  OF  H.  PS,PN  AND  PO  DEFINE  THE  LOWER  DIAGONAL,  UPPER 

C DIAGONAL  AND  MAIN  DIAGONAL  ELEMENTS  OF  V. 

C THE  COMPUTATION  IS  DONE  IN  A  DOWN-TO-UP  AND  LEFT-TO-RIGHT 

C ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  OF  COMPONENTS  OF  V  ALONG 

C THE  BOTTOM  AND  TOP  LINES  OF  THE  GRID. 

C 

DO  2  I=1,N 

C i  FOLLOWS  THE  LOWER  HORIZONTAL  LINE.  FOR  EACH  I  THE  VALUES  OF  H 

C FOR  THE  TRANSFORMED  UPPER  DIAGONAL  MATRIX  ARE  COMPUTED  ALONG  THE 

C CORRESPONDING  VERTICAL  LINE.  I  MOVES  THE  COMPUTATIONS  FROM  LEFT  TO 

C RIGHT  ALONG  THE  (BOTTOM)  HORIZONTAL  LINE.  I  IS  THE  H  INDEX  AND 

C ipN  is  THE  A2  INDEX. 

IPN=I+N 

PS=A2( I )+A2( IPN) 

PN=A2( IPN)+A2( IPN+N) 

PO=-(PS+PN) 

H(I )=-PN/(PO+W) 

C J  STEPS  UP  EACH  VERTICAL  LINE  IN  ORDER  FROM  LEFT  TO  RIGHT.  J  IS 

C THE  H  iNDEx  AND  JPN  IS  THE  A2  INDEX. 

Ll=NSO-N-N+I 
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00  1  J«IPNfLltN 

PSxPN 

JPN*J+N 

PN=A2(JPN)+A2(JPN+N) 

P0=-(PS*PN) 

H(J)=-PN/(PO+W+PS*H( J-N) ) 

1  CONTINUE 

2  CONTINUE 

C THE  VALUE  OF  H  FOR  THE  TRANSFORMED  UPPER  DIAGONAL  MATRIX 

C ALONG  THE  TOP  HORIZONTAL  LINE  OF  THE  GRID  SO  IT  NEED  NOT 

C COMPUTED  AND  STORED, 

RETURN 

END 


IS 
BE 


ZERO 


SUBROUTINE  PP2 ( N, A2 ,NSQMN, W,H,NSQ , X ,NSQP2N ) 

IMPLICIT  REAL*8( A-H,0-Z) 

DIMENSION  H(NSQMN),X(NSQ),A2(NSQP2N) 

C 

C THIS  SUBROUTINE  CONTINUES  THE  GAUSSIAN  ELIMINATION  BY  CALCULATING 

C THE  RECURSION  FORMULA  FOR  THE  RIGHT  HAND  SIDES;  I.E.  IT  CALCULATES 

C THE  RIGHT  HAND  SIDES  P(M)=X(M)  BY  THE  FOLLOWING  FORMULA: 

C P(M)=X(M)  FOR  M=1,...,N 

C P<M)=X<M)*(X(M)+PS*X(M-N) ) / ( PO-PS*H < M-N ) )  FOR  M=N+1,..,N 

C WHERE  PS,PN  AND  PO  DEFINE  THE  LOWER  DIAGONAL,  UPPER  DIAGONAL 

C AND  MAIN  DIAGONAL  ELEMENTS  OF  V  AT  POINT  M. 

C COMPUTATION  IS  DONE  IN  A  DOWN-TO-UP  AND  LEFT-TO-RIGHT  ORDER 

C oVER  THE  GRID.  FOR  EACH  VERTICAL  LINE  P(M)=X(M)  IS  CALCULATED  AT 

C THE  FIRST  POINT  AND  THEN  THE  REMAINING  POINTS  IN  ORDER  TO  TAKE 

C ADVANTAGE  OF  THE  FACT  THAT  P(M)  DOES  NOT  DEPEND  ON  A  PREVIOUS  H 

C AT  THE  FIRST  POINT  OF  ANY  VERTICAL  LINE. 

C 

DO  2  1=1, N 

C 1  FOLLOWS  THE  LOWER  HORIZONTAL  LINE.  FOR  EACH  I  THE  VALUES  OF 

C p=x  ARE  COMPUTED  ALONG  THE  CORRESPONDING  VERTICAL  LINE.  I  MOVES 

r THE  COMPUTATIONS  FROM  LEFT  TO  RIGHT  ALONG  THE  (BOTTOM)  HORIZONTAL 

C LINE.  I  IS  THE  INDEX  OF  P  (P=X)  AND  IPN  IS  THE  A2  INDEX. 

IPN=I+N 

?S=A2( I )+A2( IPN) 

PN=A2< IPN)+A2( IPN+N) 

PO=-(PS+PN> 

X(I )=X( I )/(PO+W) 

C J  STEPS  UP  EACH  VERTICAL  LINE  IN  ORDER  FROM  LEFT  TO  RIGHT.  J  IS 

c the  X  INDEX  AND  JPN  IS  THE  INDEX  OF  A2. 

L1=NSQ-N+I 

DO  1  J=IPN,L1,N 

PS  =  PN 

JPN=J+N 

PN=A2 ( JPN)+A2( JPN+N) 

PO=-(PS+PN) 

X(  J )  =  (X(  J)-PS*X(  J-N) )/(PO+W+PS*H(J-N) ) 
1  CONTINUE 
?  CONTINUE 
r.       THIS  SECTION  OE  THE  SUBROUTINE  PERFORMS  THE  BACKWARD  SUBSTITUTION 
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C FOR  THE  GAUSSIAN  ELIMINATION  PROCESS;  I.E.  IT  COMPUTES: 

C 

C X(M)»P(M)=X(M)  FOR  M*N,...,N**2-N+1 

C X(M)=P(M)+H<M)*X<M*N)=X(M)+H<M)*X<M+N)  FOR  M«N**2-N t . . . ♦ 1 

C 

I1=NSQMN+1 

DO  3  I=1,NSQMN 
C KK  IS  THE  DECREASING  INDEX  FOR  X  AND  H. 

KK=I1-I 

X(KK)=X(KK)*H(KK)*X(KK+N) 
3  CONTINUE 

RETURN 

END 


SUBROUTINE  PRODMT < N , Al , A2 ,T ,RN,NSQ ,NSQP2N ) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  AKNSQP2N)  ,A2(NSQP2N)  , T ( NSQ )  tRN(NSO) 
C 

C CALCULATE  THE  PRODUCT  M*T  WHERE  M  IS  A  MATRIX  AND  T  IS  A  VECTOR. 

C THE  COMPUTATION  IS  DONE  IN  A  LEFT-TO-RIGHT  AND  DOWN-TO-UP  ORDER 

C OVER  THE  GRID  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  OF 

C COMPONENTS  OF  M;  I.E.  THE  COMPUTATION  OF  RN=M*T  IS  PERFORMED 

C AT  THE  FIRST  GRID  POINT,  THEN  THE  INTERMEDIATE  GRID  POINTSt  AND 

C THEN  THE  LAST  GRID  POINT  OF  EACH  HORIZONTAL  LINE  IN  ORDER  TO  TAKE 

C ADVANTAGE  OF  THE  VANISHING  OF  COMPONENTS  OF  M  ALONG  THE  BOTTOM 

C AND  TOP  HORIZONTAL  LINES  AND  THE  LEFT  AND  RIGHT  VERTICAL  LINES  OF 

C THE  GRID. 

C PS,PW,PO,PE,  AND  PN  DEFINE  THE  LOWER  DIAGONAL,  LEFT  ADJACENT 

C DIAGONAL,  MAIN  DIAGONAL,  RIGHT  ADJACENT  DIAGONAL,  AND  THE  UPPER 

C DIAGONAL  OF  THE  MATRIX  M. 

C 

C DEFINE  CONSTANTS 

N1=N+1 

NN1=N1+N 

NN=N+N 

N21=NSQ-N-N+1 
C BEGIN  PRODUCT  CALCULATION 

PE=A1 (2 )+Al(3) 

PN  =  A2(N1 >+A2(NNl  ) 

PW=A1 ( I )+Al (2) 

PS  =  A2( 1 )+A2(Nl  ) 

PO=-(PE+PN+PW+PS) 

C THE  CALCULATION  OF  RN=M*T  AT  THE  FIRST  GRID  POINT  IS  PERFORMED 

C HERE.  NOTE  THE  COMPONENT  PS  OF  M  VANISHES  ALONG  THE  BOTTOM 

C HORIZONTAL  LINE. 

RN(1 )=PO*T(l )+PE*T(2)+PN*T(Nl) 

DO  15  1=2, N 

C 1  MOVES  THE  COMPUTATION  ALONG  THE  REMAINING  GRID  POINTS  ON  THE 

C BOTTOM  HORIZONTAL  LINE.  RN=M*T  AT  THESE  POINTS  IS  COMPUTED  HERE.  I 

C IS  THE  BASIS  FOR  THE  INDICES  OF  Al,  A2,  AND  RN. 

PE  =  A1  (  I  +  D+AK  1  +  2) 

PN=A2( I+N)+A2( I+NN) 

PW  =  A1  (  I  )+AK  1  +  1) 


72 
PS=A2(I )+A2< I+N) 
PO=-(PE+PN+PW+PS) 

C F0R  I=N  THE  VALUE  OF  RN(N)  IS  INCORRECT  AND  IS  CORRECTED  BELOW, 

RNU  )=PW*T(  I-l)+PO*T(  I  )+PE*T(I+l)+PN*T(I+N) 
15  CONTINUE 

RN=M*T  AT  THE  LAST  GRID  POINT  ON  THE  BOTTOM  HORIZONTAL  LINE  IS 

CALCULATED  HERE  SEPARATELY  FROM  THE  INTERMEDIATE  POINTS  ON  THE 

LINE  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  COMPONENT  PE  OF 

M.  THIS  IS  DONE  BY  SUBTRACTING  PE*T<N+1>  FROM  THE  RN(N)  JUST 

CALCULATED. 

RN(N)=RN(N)-PE*T(N+1  ) 
N1=N+1,  N21=N**2-2*N+1 

M1=N1 

M2  =  N 

DO  30  J=N1,N21,N 

J  MOVES  THE  COMPUTATION  UP  TO  THE  NEXT  HORIZONTAL  LINE  OF  THE 

GRID.  FOR  EACH  J  RN=M*T  IS  CALCULATED  AT  THE  FIRST  GRID  POINT,  THE 

INTERMEDIATE  POINTS,  AND  THE  LAST  GRID  POINT  OF  THE  GIVEN 

HORIZONTAL  LINE.  Ml,  M2  AND  J  ARE  THE  BASES  FOR  THE  INDICES  OF  Alt 

A2,  AND  RN.  RN=M*T  FOR  THE  FIRST  GRID  POINT  OF  THE  HORIZONTAL  LINE 

CALCULATED  HERE  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING 

COMPONENT  PW  OF  M. 

Ml=Ml+3 

M2=M2+1 

PE  =  A1 (Ml )+Al (Ml  +  1  ) 

PN=A2(M2 )+A2(M2+N) 

PW=A1 (Ml-1 )+Al (Ml) 

PS=A2(M2-N)+A2(M2) 

PO=-(PE+PN+PW+PS) 

RN( J)=PS*T( J-N)+PO*T( J)+PE*T( J  +  l )+PN*T( J+N) 

J1=J+1 

J2=J+N-1 

DO  25  I=J1,J2 

1  MOVES  THE  COMPUTATION  ALONG  THE  REMAINING  GRID  POINTS  ON  THE 

GIVEN  HORIZONTAL  LINE  AND  RN=M*T  FOR  THOSE  POINTS  IS  COMPUTED 

HERE.  Ml,  M2  AND  I  ARE  THE  BASES  FOR  THE  INDICES  OF  Al,  A2  AND  RN. 

M1=M1+1 

M2=M2+1 

PE  =  A1  (Ml  )+AKMl  +  l) 

PN=A2(M2)+A2(M2+N) 

PW=A1 (Ml-1 )+Al (Ml) 

PS=A2(M2-N)+A2(M2) 

PO=-(PE+PN+PW+PS) 

RN(  I  )=PS*T(  I-N)+PW*T( I-1)+P0*T( I )+PE*T( 1+1  )+PN*T(  I+N) 
25  CONTINUE 

THE  COMPUTATION  OF  RN=M*T  AT  THE  LAST  GRID  POINT  ON  THE  GIVEN 

HORIZONTAL  LINE  IS  DONE  HERE  SEPARATELY  FROM  THE  OTHER  GRID  POINTS 

ON  THE  LINE  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  COMPONENT 

PE  OF  M.  THIS  IS  DONE  BY  SUBTRACTING  PE*T(J2+l)  FROM  THE  RN(J2) 

JUST  CALCULATED. 

PN(J2)=RN(J2)-PE*T(J2+1) 
30  CONTINUE 

L=J2+1 

L  MOVES  THE  COMPUTATION  UP  TO  THE  FIRST  GRID  POINT  ON  THE  TOP 

HORIZONTAL  LINE  AND  RN=M*T  AT  THAT  POINT  IS  COMPUTED  HERE.  Ml,  M2 

AND  L  ARE  THE  BASES  FOR  THE  INDICES  OF  Al,  A2  AND  RN.  NOTE  THE 

COMPONENT  PN  OF  M  VANISHES  ALONG  THIS  LINE  AND  PW  VANISHES  AT 
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THIS    POINT. 

Ml=Ml+3 

M2=M2+1 

PE  =  A1  (Ml  l+AKMl  +  1) 

PN=A2(M2)+A2(M2+N) 

PW=A1 (Ml-1 I+A1IM1) 

PS=A2(M2-N)+A2(M2) 

PO=-(PE+PN+PS*PW) 

RN(L)=PS*T(L-N)+PO*T(L)+PE*T<L+l) 

L1=L+1 

DO  45  K=L1»NSQ 

K  MOVES  THE  COMPUTATION  ALONG  THE  REMAINING  GRID  POINTS  ON  THE 

TOP  HORIZONTAL  LINE  AND  RN=M*T  AT  THOSE  POINTS  IS  COMPUTED  HERE. 

Ml,  M2  AND  K  ARE  THE  BASES  FOR  THE  INDICES  OF  Alt  A2  AND  RN. 

M1=M1+1 

M2=M2+1 

PE  =  A1  (Ml  )+AKMl+l) 

PN=A2(M2 )+A2(M2+N) 

PW=A1(M1-1  )+AKMl) 

PS=A2(M2-N)+A2(M2) 

PO=-(PE+PN+PW+PS) 
FOR  K  =  NSQ  THE  VALUE  OF  RN(NSQ)  IS  COMPUTED  BELOW. 

IF(K.EQ.NSO)  GO  TO  40 

RN(K)=PS*T(K-N)+PW*T(K-1)+P0*T(K)+PE*T(K+1) 

GO  TO  45 

THE  COMPUTATION  RN=M*T  AT  THE  LAST  GRID  POINT  ON  THE  TOP 

HORIZONTAL  LINE  IS  PERFORMED  HERE  SEPARATELY  FROM  THE  OTHER  GRID 

POINTS  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  OF  THE 

COMPONENTS  PE  AND  PN  OF  M. 

40  RN(K)=PS*T(K-N)+PW*T(K-1 )+PO*T(K) 
45  CONTINUE 

RETURN 

END 


SUBROUTINE  RES IDL ( RN »0 »NSO ) 
IMPLICIT  REAL*8(A-H,0-Z ) 
DIMENSION  RN(NSO) »0(NSQ) 

•THIS  SUBROUTINE  CALCULATES  THE  RESIDUAL  RN=Q-M*T=Q-RN  WHERE  0  IS 
AN  INPUT  VECTOR  AND  RN=M*T  IS  A  VECTOR  COMPUTED  BY  THE  SUBROUTINE 
■PRODMT. 

DO  1  1=1, NSQ 

RN( I )=0(I )-RN( I ) 

CONTINUE 

RETURN 

END 
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APPENDIX  C 
NEUMANN  SUBROUTINES  FOR  THE  STONE  AND  ADI  PROGRAMS 
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SUBROUTINE  ALTERM ( N, Al , A2, BLtCL t DL , EU, FU, ALPHA, NSQ,NSQP2N, NSOMN ) 

IMPLICIT  REAL*8( A-H,0-Z) 

DIMENSION  A1(NSQP2N),A2(NSQP2N) , BL ( NSQMN ) ,CL < NSQMN ) , DL ( NSQ ) , EU < NSQ 
XMN) tFU(NSQMN) 
C 

C CALCULATE  THE  ALTERED  MATRIX  M+N=L*U  FROM  THE  GIVEN  MATRIX  M.  L  IS 

C A  TRI-DIAGONAL  MATRIX  CONSISTING  OF  THE  MAIN  DIAGONAL,  THE  LEFT 

C ADJACENT  DIAGONAL,  AND  A  LOWER  DIAGONAL  AT  A  DISTANCE  N  FROM  THE 

C MAIN  ONE.  THE  ELEMENTS  OF  THESE  DIAGONALS  ARE  STORED  IN  THE 

C VECTORS  DL,  CL,  AND  BL  RESPECTIVELY.  U  IS  A  TRI-DIAGONAL  MATRIX 

C CONSISTING  OF  THE  MAIN  DIAGONAL  EVERY  ELEMENT  OF  WHICH  IS  EQUAL 

C TO  THE  CONSTANT  ONE,  THE  RIGHT  ADJACENT  DIAGONAL,  AND  AN  UPPER 

c DIAGONAL  AT  A  DISTANCE  N  FROM  THE  MAIN  ONE,  THE  ELEMENTS  OF  THE 

C THE  TWO  OFF  DIAGONALS  ARE  STORED  IN  THE  VECTORS  EU  AND  FU 

C RESPECTIVELY.  VANISHING  ELEMENTS  OF  THE  DIAGONALS  ARE  NOT  STORED 

C jN  THESE  AFORE  MENTIONED  VECTORS. 

C COMPUTATION  IS  DONE  IN  A  LEFT-TO-RIGHT  AND  DOWN-TO-UP  ORDER 

C SINCE  THE  ELEMENTS  OF  L*U  CALCULATED  AT  ANY  GRID  POINT  ONLY 

C DEPEND  ON  PREVIOUSLY  CALCULATED  ELEMENTS  AT  PREVIOUS  GRID  POINTS. 

c PS,PW,PO,PE,  AND  PN  DEFINE  THE  LOWER  DIAGONAL,  LEFT  ADJACENT 

C DIAGONAL,  MAIN  DIAGONAL,  RIGHT  ADJACENT  DIAGONAL,  AND  UPPER 

C DIAGONAL  OF  THE  MATRIX  M. 

C 

C DEFINE  CONSTANTS 

N1=N+1 

NN1=N1+N 

NN=N+N 

N21=NSQ-N-N+1 

NSOMl=NSO-l 
C GET  M+N=L*U  FROM  M 

PE=A1(2)+A1(3) 

PN=A2(N1 )+A2(NNl ) 

PW=A1 ( 1 )+Al<2) 

PS=A2(1)+A2(N1) 

PO=-(PE+PN+PW+PS) 

C THE  COMPUTATION  IS  PERFORMED  AT  THE  FIRST  POINT,  THEN  THE 

C INTERMEDIATE  POINTS,  AND  THEN  THE  LAST  POINT  OF  THE  GIVEN 

C HORIZONTAL  LINE  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  OF 

C COMPONENTS  OF  M  ALONG  THE  BOTTOM  AND  TOP  HORIZONTAL  LINES,  AND 

C THE  LEFT  AND  RIGHT  VERTICAL  LINES  OF  THE  GRID.  THE  ELEMENTS  OF  L*U 

C FOR  THE  FIRST  GRID  POINT  ARE  COMPUTED  HERE. 

DL(1 )=PO 
C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

EU(1 )=(PE+PE)/DL( 1) 

FU(1 )=<PN+PN)/DL(1) 

DO  15  1=2, N 

C 1  FOLLOWS  THE  REMAINING  GRID  POINTS  ON  THE  BOTTOM  HORIZONTAL 

C LINE.  THE  ELEMENTS  OF  L*U  FOR  THOSE  POINTS  ARE  COMPUTED  HERE.  I  IS 

C THE  BASIS  FOR  THE  INDICES  OF  Al,  A2,  AND  THE  ELEMENTS  OF  L*U. 

PE  =  A1  (I  +  D+AK  1  +  2) 

PN=A2( I+N)+A2( I+NN) 

PW=A1  (  I  )+AK  1  +  1) 

PS  =  A2( I  )+A2(I+N) 

PO=-( PE+PN+PW+PS) 

CL( 1-1 )=PW/(1.0DOO+ALPHA*FU(I-l ) ) 

DL( I )=PO+CL( I-l)*( ALPHA*FU( I-1)-EU( 1-1 ) ) 
C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 
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FUU  )=(PN+PN-ALPHA*CL(I-1)*FU<I-1)  )/DL(I) 

C. FOR  THE  LAST  GRID  POINT  ON  THE  BOTTOM  HORIZONTAL  LINE  THE  RIGHT 

C ADJACENT  DIAGONAL  ELEMENT  OF  L*U  VANISHES  SO  THE  COMPUTATION  OF 

C EU(N)  IS  INCORRECT  AND  WILL  BE  CORRECTED  IN  STATEMENT  16. 

EU(I )=PE/DL(I ) 

15  CONTINUE 

16  EU(N)=O.ODOO 

C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

CL(N-l)=(PW+PW)/(1.0DOO+ALPHA*FU(N-l)) 

C L  AND  M  ARE  COUNTER  VARIABLES  THAT  AID  IN  CALCULATING  THE  CORRECT 

C SUBSCRIPTS  FOR  THE  COMPONENTS  OF  L*U.  M2  INCREMENTS  THE  Al  INDEX 

c AT  EACH  HORIZONTAL  LINE.  N1*N*1 1 N21»N**2-2*N+1  . 

M=0 

M2  =  0 

L  =  2 

DO  30  J=N1,N21,N 

c j  MOVES  THE  COMPUTATION  UP  TO  THE  NEXT  HORIZONTAL  LINE  OF  THE 

c GRID.  FOR  EACH  J  CALCULATIONS  ARE  DONE  AT  THE  FIRST  POINT,  THE 

c INTERMEDIATE  POINTS,  AND  THE  LAST  POINT  OF  THE  GIVEN  HORIZONTAL 

c LINE.  J  IS  THE  BASIS  FOR  THE  INDICES  OF  Al,  A2,  AND  THE  ELEMENTS 

c 0F  L*u.  THE  CALCULATIONS  FOR  THE  FIRST  POINT  ON  EACH  HORIZONTAL 

C LINE  ARE  DONE  HERE. 

M2=M2+2 

PE=A1( J+M2+1 )+AK  J  +  M2  +  2) 

PN=A2(J+N)+A2< J+NN) 

PW=A1( J+M2 )  +  Al( J+M2+1) 

PS=A2(J)+A2( J  +  N) 

PO=-(PE+PN+PW+PS> 

BL(J-N)=PS/( 1.0D00+ALPHA*EU< J-N-M) ) 

DL(J)=PO+BL( J-N)*(ALPHA*EU( J-N-M ) -FU ( J-N ) ) 
C****»THIS  CHANGE  IS  REOUIRED  FOR  THE  NEUMANN  PROBLEM. 

EU(J-1-M)=(PE+PE-ALPHA*BL( J-N )*EU< J-N-M) )/DL(J) 

FU(J)=PN/DL( J) 

J1=J+1 

J2=J+N-2 

DO  25  I=J1,J2 

c i  MOVES  THE  COMPUTATION  ALONG  THE  REMAINING  GRID  POINTS  ON  THE 

C GIVEN  HORIZONTAL  LINE  AND  THOSE  COMPUTATIONS  OF  THE  ELEMENTS  OF 

c L*U  are  DONE  HERE.  I  IS  THE  BASIS  FOR  THE  INDICES  OF  Al,  A2,  AND 

c the  ELEMENTS  OF  L*U. 

PE  =  A1  (  I+M2  +  D+AK  I+M2+2) 

PN=A2( I+N)+A2( I+NN) 

PW=A1 (I+M2)+A1( I+M2+1) 

PS  =  A2(  I  )+A2( I+N) 

PO=-(PE+PN+PW+PS) 

BL( I-N)=PS/( 1.0D00+ALPHA*EU( I-N-M) ) 

CL( I-L)=PW/( 1.0D00+ALPHA*FU( 1-1) ) 

DL(I )=P0+BL(I-N)*(ALPHA*EU(I-N-M)-FU(I-N))+CL(I-L)*(ALPHA*FU(I-1) 

XEU( I-2-M) ) 
EU(I-1-M)=(PE-ALPHA*BL( I-N)*EU( I-N-M) )/DL(I ) 
FU( I )=(PN-ALPHA*CL( I-L)*FU( 1-1) )/DL( I ) 
25  CONTINUE 
I=.J2  +  1 

C THE  COMPUTATION  OF  THE  ELEMENTS  OF  L*U  FOR  THE  LAST  POINT  ON  THE 

c GIVEN  HORIZONTAL  LINE  ARE  DONE  HERE  SEPARATELY  FROM  THE  OTHER  GRID 

C POINTS  ON  THE  LINE  IN  OROER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING 

C OF  THE  COMPONENT  EU  OF  L*U  AT  THIS  POINT. 
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BL(J-1)=PS 
C*«#**THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 
CL(I-L)  =  (PW  +  PW)/(1.0DOO-»-ALPHA*FU(I-l)) 

DLU >=PO+BL< I-N)*(-FU< I-N) ) +CL ( I-L >*( ALPHA*FU( 1-1 ) -EU ( I-2-M ) ) 
FU( I )=<PN-ALPHA*CL( I-L)*FU(I-1) )/DL(I ) 
M  =  M  +  1 
L  =  L  +  1 
30  CONTINUE 

C 1  MOVES  THE  COMPUTATION  UP  TO  THE  FIRST  POINT  ON  THE  TOP 

c HORIZONTAL  LINE  AND  THE  ELEMENTS  OF  L*U  FOR  THAT  POINT  ARE 

C CALCULATED  HERE.  I  FORMS  THE  BASIS  FOR  THE  INDICES  OF  Alt  A2»  AND 

C THE  ELEMENTS  OF  L*U. 

PE=A1(I+M2+1)+A1( I+M2+2) 

PN  =  A2U+N)*A2(  I+NN) 

PW=A1( I+M2)+A1( I+M2+1) 

PS=A2< I J+A2U+N) 

PO=-(PE+PN+PW+PS) 
C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

BL(I-N)  =  (PS  +  PS)/(  1.0D00+ALPHA*EUU-N-M)  ) 

DL( I )=PO+BL( I-N)*( ALPHA*EU( I-N-M )-FU ( I-N ) ) 
C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

EU( I-1-M)=(PE+PE-ALPHA*BL< I-N ) *EU ( I-N-M » )/DL( I ) 

11=1+1 

DO  45  K=I1,NSQM1 

C K  MOVES  THE  COMPUTATION  ALONG  THE  REMAINING  GRID  POINTS  ON  THE 

C TOP  HORIZONTAL  LINE  AND  THOSE  COMPUTATIONS  OF  THE  ELEMENTS  OF  L*U 

C ARE  DONE  HERE.  K  FORMS  THE  BASIS  FOR  THE  INDICES  OF  Al,  A2 »  AND 

C THE  ELEMENTS  OF  L*U. 

PE=A1 (K+M2+1 )+Al(K+M2+2) 

PN=A2(K+N)+A2(K+NN) 

PW=A1 (K+M2)+A1(K+M2+1) 

PS=A2(K )+A2(K+N) 

PO=-(PE+PN+PW+PS) 
C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

BL(K-N)=(PS+PS)/( 1.0D00+ALPHA*EU(K-N-M) ) 

CL(K-N)=PW 

DL(K)=PO+BL(K-N)*( ALPHA*EU ( K-N-M ) -FU( K-N ) )+CL < K-N ) *(-EU ( K-N ) ) 

EU(K-N+1 )=(PE-ALPHA*BL( K-N) *EU( K-N-M) )/DL(K) 
45  CONTINUE 

C THE  COMPUTATION  OF  THE  ELEMENTS  OF  L*U  FOR  THE  LAST  GRID  POINT 

C OF  THE  TOP  HORIZONTAL  LINE  ARE  DONE  HERE  SEPARATELY  FROM  THE  OTHER 

C GRID  POINTS  ON  THE  LINE  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE 

C VANISHING  OF  THE  COMPONENTS  EU  AND  FU  OF  L*U  AT  THIS  POINT. 

C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

BL(NSQMN)=PS+PS 

CL (NSQMN)=PW+PW 

DL(NSQ)=PO+BL(NSQMN)*(-FU(NSQMN) ) +CL ( NSQMN )* ( -EU ( NSQMN ) ) 

RETURN 

END 
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SUBROUTINE  BALTRM(N, Al, A2, B,C,D, E,F, ALPHA, NSQ, NSOP2N, NSQMM ) 
IMPLICIT  REAL*8( A-H,0-Z) 

01  HENS  I  ON  AKNSQP2N),A2<NSQP2N),B<NSQMN),C<NSQ),D(NSQUE<NSQ),F<N 

XQMN) 

C 

C CALCULATE  THE  ALTERED  HATRIX  M+N«L*U  FROH  THE  GIVEN  MATRIX  M.  L  IS 

C A  TRI-OIAGONAL  MATRIX  CONSISTING  OF  THE  MAIN  DIAGONAL,  THE  LEFT 

C ADJACENT  DIAGONAL,  AND  A  LOWER  DIAGONAL  AT  A  DISTANCE  N  FROM  THE 

C MAIN  ONE.  THE  ELEMENTS  OF  THESE  DIAGONALS  ARE  STORED  IN  THE 

C VECTORS  D,  C,  AND  B  RESPECTFULLY.  U  IS  A  TRI-DIAGONAL  MATRIX 

C CONSISTING  OF  THE  MAIN  DIAGONAL  EVERY  ELEMENT  OF  WHICH  IS  EQUAL 

C TO  THE  CONSTANT  ONE,  THE  RIGHT  ADJACENT  DIAGONAL,  AND  AN  UPPER 

C DIAGONAL  AT  A  DISTANCE  N  FROM  THE  MAIN  ONE.  THE  ELEMENTS  OF  THE 

C THE  TWO  OFF  DIAGONALS  ARE  STORED  IN  THE  VECTORS  E  AND  F 

C RESPECTFULLY.  VANISHING  ELEMENTS  OF  C  AND  E  ARE  STORED  AS  ZEROES 

C IN  THESE  AFORE  MENTIONED  VECTORS. 

C THE  COMPUTATION  IS  PERFORMED  IN  REVERSE  ORDER  FROM  SUBROUTINE 

C ALTERM.  THAT  IS,  THE  SOLUTION  IS  COMPUTED  AT  THE  GRIO  POINTS 

C IN  THE  J  =  N  ROW  WITH  I  ASSUMING  THE  VALUES  1  UP  TO  N.  THEN 

C j  is  DECREMENTED  BY  1  AND  I  AGAIN  ASSUMES  THE  RANGE  OF 

C VALUES  1  THROUGH  N.  THIS  LEFT-TO-RIGHT  AND  UP-TO-DOWN 

C ORDERING  OF  THE  GRID  POINTS  FOR  EVEN  ITERATIONS  BEGINNING  WITH 

C THE  SECOND  MAY  BE  CONSIDERED  A  REORDERING  OF  THE  ORIGINAL  POINTS. 

C PS,PW,PO,PE,  AND  PN  DEFINE  THE  LOWER  DIAGONAL,  LEFT  ADJACENT 

C DIAGONAL,  MAIN  DIAGONAL,  RIGHT  ADJACENT  DIAGONAL,  AND  UPPER 

C DIAGONAL  OF  THE  MATRIX  M. 

C 

C COMPUTE  CONSTANTS 

NM2=N-2 

NP1=N+1 

NSQP1=NSQ+1 

NSOPN=NSO+N 

NM1SQ=NSQ-N-N+1 

NSQMN1=NSQMN+1 
C THE  SUBSCRIPTS  OF  Al  AND  A2  FOLLOW  THE  NORMAL  ORDERING. 

PE=A1 (NSQPN)+A1(NSQPN+1) 

PN=A2(NSQP1)+A2(NSQP1+N) 

PW  =  A1  (NSOPN-D+AKNSOPN) 

PS=A2(NSQP1-N)+A2(NSQP1) 

PO=-(PE+PN+PW+PS) 

D(l )=P0 
C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

E(l )=(PE+PE)/D(1) 

F(1)=(PS+PS)/D(1) 

C(l )=0.0D00 

C OBSERVE  THAT  THE  MATRIX  M  IS  REORDERED  WITH  THE  RESULT  THAT  FOR 

C A  GIVEN  ROW  I,  PN  AND  PS  ARE  THE  COEFFICIENTS  OF  X(I-N)  AND  X  (  I+n) 

C RESPECTFULLY.  THIS  REVERSES  THEIR  ROLE  IN  THE  OLD  ORDERING  WHERE 

c PN  AND  ps  ARE  THE  COEFFICIENTS  OF  XU+N)  AND  XU-N). 

C THE  ROLE  OF  PE  AND  PW  IS  UNCHANGED. 

C THE  FOLLOWING  LOOP  COMPUTES  C,D,E,  AND  F  ALONG  THE  TOP 

C HORIZONTAL  LINE  OF  THE  GRID. 

DO  10  1=2, N 

C Lfjnp  PARAMETERS  FOLLOW  GRID  POINTS  IN  THE  NEW  ORDERING.  THESE 

C PARAMETERS  AVOID  SUPERFLUOUS  ADDITIONS  OF  I. 

1*1*1-1 

NSQI-NSQ-MM1 
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NSQP1I*NSQP1+IM1 

NSQPNI»NSQPN*IM1 

PE»A1(NSQPNI )  +  AKNSQPNI  +  l) 

PN=A2(NSQP1I )+A2<NSQPlI+N) 

PW=A1(NSQPNI-1)+A1(NSQPNI) 

PS»A2(NSQPlI-N)«-A2<NSQPin 

PO=-(PE+PN*PW+PS) 

C(I  )«PW/(1.0D00+ALPHA*F<  I  — 1 1  > 

D(I  )«P0+C(I)*(ALPHA*FU-1)-EU-1>) 

E(I )»PE/D(  I  ) 

C FOR  I»N  THIS  VALUE  OF  E  IS  INCORRECT  AND  IS  CORRECTED  IN 

C STATEMENT  11. 

C»****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

F(I  )  =  (PS  +  PS-ALPHA*C( I)*F< I-1))/D(I ) 

10  CONTINUE 

11  E(N)=O.ODOO 

C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

C(N)=(PW+PW)/<1.0D00+ALPHA*F(N-1) ) 

C THE  NEXT  SEGMENT  COMPUTES  B,C,D,E,  AND  F  FROM  THE  SECOND 

C HORIZONTAL  LINE  (IN  UP-TO-DOWN  ORDER)  TO  AND  INCLUDING  THE  N-l 

C HORIZONTAL  LINE.  THE  SEGMENT  IS  A  NESTED  LOOP.  ONE  EXECUTION  OF 

C THE  OUTER  LOOP  COMPUTES  THESE  COMPONENTS  ALONG  THE  FIXED 

C HORIZONTAL  LINE.  CERTAIN  PARAMETERS  AND  COMPONENTS  ARE  ZERO  AT  THE 

C ENDPOINTS  OF  HORIZONTAL  LINES.  ACCORDINGLY,  THE  COMPUTATIONS 

C FOR  A  FIXED  HORIZONTAL  LINE  PROCEED  AS  FOLLOWS: 

C (1)  C  =  0  AND  PW=0  AT  THE  LEFT  ENDPOINT.  THIS  FACT  ABOUT  PW 

C REQUIRES  THE  COMPUTATION  AT  THIS  POINT  TO  BE  A  SPECIAL  CASE 

C N0TE  THAT  SINCE  C*0  AT  THE  ENDPOINTS,  IT  IS  POSSIBLE  TO 

C AVOID  RESERVING  CORRESPONDING  STORAGE  IN  THE  C  ARRAY.  THIS 

C WAS  NOT  DONE  IN  ORDER  TO  SIMPLIFY  THE  CALCULATION  OF 

C SUBSCRIPTS  OF  C. 

C (2)  POINTS  2  THROUGH  N-l  ON  THE  GIVEN  HORIZONTAL  LINE  ARE 

C COMPUTED  ITERATIVELY. 

c (3)  E  =  0  AND  PE  =  0  AT  THE  right  ENDPOINTS.  PE  =  0  MAKES  THIS  A 

C SPECIAL  CASE.  NOTE  ALSO,  THAT  E  =  0  AT  THESE  ENDPOINTS  MAKES 

C IT  POSSIBLE  TO  AVOID  RESERVING  CORRESPONDING  STORAGE  IN  THE 

C E  ARRAY.  THIS  WAS  NOT  DONE  IN  ORDER  TO  SIMPLIFY  THE 

C CALCULATION  OF  THE  SUBSCRIPTS  OF  E. 

C (4)  NOTE  THAT  THE  VALUES  OF  C  AND  E  THAT  ARE  ZERO  AT  CERTAIN 

C ENDPOINTS  ARE  STORED  AS  ZERO  SO  THAT  IF  ONE  DESIRES  A 

C PRINTOUT  OF  THE  ARRAYS  FOR  DEBUGGING  PURPOSES  THE  CORRECT 

C VALUES  OF  C  AND  E  WILL  BE  STORED  FOR  THE  CORRESPONDING 

C ENDPOINTS. 

C THE  OUTER  LOOP  MOVES  THE  COMPUTATION  FROM  ONE  HORIZONTAL  LINE  TO 

C THE  NEXT  LOWER. 

Kl=-N-2 
C Kl  STEPS  Al  DOWN 

K2=-N 

C K2  STEPS  A2  DOWN.  THE  LOOP  VARIABLE  FOLLOWS  THE  INDEX  OF  POINTS 

C OF  THE  LEFT  VERTICAL  LINE  IN  THE  NEW  ORDERING. 

DO  30  I=NP1,NM1SQ,N 

Kl=KH-N+2 

K2=K2+N 

Ll=NSQ-2-Kl 

C L1  is  THE  INDEX  OF  Al  AT  THE  LEFT  GRID  POINT  OF  THE  CURRENT 

C HORIZONTAL  LINE. 

L2=NSQ-N+1-K2 


80 

•L2  IS  THE  INDEX  OF  A2  AT  THE  LEFT  GRID  POINT  OF  THE  CURRENT 

-HORIZONTAL  LINE. 
PE»A1<L1)*AKL1  +  1) 
PN=A2(L2H-A2(L2-*-N> 
PW«Al<Ll-lH-Ai(Ll> 

PS»A2(L2-N)+A2(L2) 

PO=-(PE*PN+PW«-PS) 

IMN»I-N 

B(IMN)«PN/<1.0DOO+ALPHA*EUMN>) 

C(I)«0.0D00 

D(I )«PO*B( IMN)*(ALPHA*E(IMN)-F(IMN>) 
C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

E(I)=<PE+PE-ALPHA*B(IMN)*E(IMN) >/D<I) 

FU)=PS/D(I) 

C jj  FOLLOWS  THE  INDEX  OF  THE  INTERIOR  GRID  POINTS  IN  THE  NEW 

C ORDERING. 

DO  20  J  =  UNM2 

JJ=I+J 

L1J«L1+J 

L2J=L2+J 

PE=A1(L1J)+A1<L1J+1) 

PN=A2(L2J)+A2(L2J+N) 

PW=AKL1J-1>+A1(L1J> 

PS=A2(L2J-N)+A2(L2J) 

PO=-(PE+PN+PW+PS) 

JMN=JJ-N 

BUMN)=PN/(1.0+ALPHA*E<  JMN) ) 

C(JJ)=PW/(1.0D00+ALPHA*F( J  J-l  > ) 

D(JJ)  =  PO+B<JMN>*<ALPHA*E(JMN)-F<JMN>H-C(JJ)*ULPHA*F<JJ-l)-EUJ-l) 

X) 
E(JJ)=(PE-ALPHA*B< JMN)*E( JMN) )/D(JJ) 
F(JJ)=(PS-ALPHA*C( JJ)*F( JJ-1))/D(JJ) 
20  CONTINUE 
c this  SECTION  COMPUTES  THE  COMPONENTS  AT  THE  RIGHT  ENDPOINT  ON  THE 

C GIVEN  HORIZONTAL  LINE. 

L1J=L1J+1 

L2J=L2J+1 

PE  =  Al(LU)-»-Al(LlJ+l) 

PN=A2(L2J)*A2<L2J+N) 

PW  =  A1  (LIJ-D+AKLIJ) 

PS=A2(L2J-N)+A2(L2J) 

PO=-(PE+PN+PW+PS) 
c E=0  AND  PE=0  AT  THIS  ENDPOINT.  ZERO  IS  STORED  IN  E(JJ). 

JJ=JJ+1 

JMN=JJ-N 

B( JMN)=PN 

C SINCE  E  =  0 

C*»»*»THIS  CHANGE  IS  REOUIRED  FOR  THE  NEUMANN  PROBLEM. 

C(JJ)=(PW+PW)/( 1.0D00+ALPHA*F(JJ-1) ) 

D(JJ)=PO+B(JMN)*(-F(JMN))+C( J J ) * ( ALPHA*F ( J J-l ) -E ( J J-l > ) 

E(JJ  )  =  0.0D00 
C SINCE  PW  DOES  NOT  APPEAR. 

FUJ  )=(PS-ALPHA*C< JJ)*F< JJ-1)  )/D<  J  J  ) 
30  CONTINUE 
C THE  COMPUTATION  IS  NOW  AT  THE  FINAL  HORIZONTAL  LINE. 

PF  =  A1 (2  H-Al (3) 

PN=A2(NPl H-A2(NP1+N) 
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PW«A1U)«-A1<2) 

PS»A2(NP1-N)«-A2<NP1) 

PO*-(PE*PN+PW+PS) 

C NM1SQ=(N**2)-(2*N)*1 

C**«»*THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

fi(NMlSQ)s<PN+PN)/U.ODOO+ALPHA*E<NMlSQH 
c NSQMN 1«  <  N**2  )  -N*  1 

D<NSQMN1  )*P0+B(NM1SQ)*(  ALPMA*E<  NM1SO-F  (NM1SQ)  ) 
C*****THIS    CHANGE    IS    REQUIRED    FOR    THE    NEUMANN    PROBLEM. 

E(NSQMN1)«(PE*PE-ALPHA*B<NM1SQ)*E(NM1SQ) )/D<NSQMNl> 

C(NSQMN1)=0.0D00 

C THE  COMPUTATION  AT  THE  INTERIOR  POINTS  OF  THE  FINAL  HORIZONTAL 

C LINE  IS  DONE  HERE.  THE  LOOP  VARIABLE  INDEXES  THE  GRID  POINTS 

c im  The  NEW  ORDERING.  K  ADJUSTS  THE  INDICES  Kl  AND  K2  OF  Al  AND  A2 

K  =  -l 

L1=NSQMN1+1 

L2=NSQ-1 

DO    60    J=LltL2 

K=K  +  1 

Kl=3+K 

K2=N+2+K 

PE=A1(K1)+A1(K1+1> 

PN=A2<K2)+A2(K2+N) 

PW  =  AKK1-1 )+Al(Kl) 

PS=A2(K2-N)+A2(K2) 

PO=-(PE+PN+PW+PS) 

JMN=J-N 
C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

B(JMN)=(PN+PN)/(  1.0D00+ALPHA*E( JMN)  ) 

C(J)=PW 

D(J)=PO+B( JMN)*( ALPHA*E( JMN)-F( JMN) ) +C ( J >*(-E( J-l ) ) 

E  (J  )=(PE-ALPHA*B( JMN)*E( JMN) )/D( J) 
60  CONTINUE 
C WE  ARE  NOW  AT  THE  RIGHT  ENDPOINT.  NSQMN= ( N**2 )-N. 

K1=N  +  1 

K2=N+N 

PE  =  A1  (Kl  )+AL(Kl+l) 

PN=A2(K2)+A2(K2+N) 

PW  =  A1  (K1-1I+AHK1) 

PS=A2(K2-N)+A2(K2) 

PO=-(PE+PN+PW+PS) 
C***«*THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

B(NSQMN)=PN+PN 

C<NSQ)=PW+PW 

D(NSQ)=PO+B(NSQMN)*(-F (NSQMN) ) +C ( NSQ ) * < -E ( NSQ-1 ) ) 

E(NSQ)=O.ODOO 

RETURN 

END 
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SUBROUTINE  PRODMT< N, Alt A2, T,RN,NSQ,NSQP2N ) 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A1(NSQP2N),A2<NSQP2N),T(NSQ),RN(NSQ) 
C 

C CALCULATE  THE  PRODUCT  M*T  WHERE  M  IS  A  MATRIX  AND  T  IS  A  VECTOR. 

C THE  COMPUTATION  IS  DONE  IN  A  LEFT-TO-RIGHT  AND  DOWN-TO-UP  ORDER 

C OVER  THE  GRID  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  OF 

C COMPONENTS  OF  M;  I.E.  THE  COMPUTATION  OF  RN«M*T  IS  PERFORMED 

C AT  THE  FIRST  GRID  POINT,  THEN  THE  INTERMEDIATE  GRID  POINTSt  AND 

C THEN  THE  LAST  GRID  POINT  OF  EACH  HORIZONTAL  LINE  IN  ORDER  TO  TAKE 

C ADVANTAGE  OF  THE  VANISHING  OF  COMPONENTS  OF  M  ALONG  THE  BOTTOM 

C AND  TOP  HORIZONTAL  LINES  AND  THE  LEFT  AND  RIGHT  VERTICAL  LINES  OF 

C THE  GRID. 

C PS,PW,PO,PE,  AND  PN  DEFINE  THE  LOWER  DIAGONAL*  LEFT  ADJACENT 

C DIAGONAL,  MAIN  DIAGONAL,  RIGHT  ADJACENT  DIAGONAL,  AND  THE  UPPER 

C DIAGONAL  OF  THE  MATRIX  M. 

C 

C DEFINE  CONSTANTS 

N1*N+1 

NN1=N1+N 

NN-N-t-N 

N21=NSQ-N-N+1 
C BEGIN  PRODUCT  CALCULATION 

PE  =  AK2)+A1(3) 

PN=A2(N1 )+A2(NNl) 

PW=A1(1)+A1<2) 

PS=A2(1)+A2(N1) 

PO=-(PE+PN+PW+PS) 

C THE  CALCULATION  OF  RN=M*T  AT  THE  FIRST  GRID  POINT  IS  PERFORMED 

C HERE.  NOTE  THE  COMPONENT  PS  OF  M  VANISHES  ALONG  THE  BOTTOM 

C HORIZONTAL  LINE. 

C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

RN(1 )=PO*T(l)+(PE+PE)*T(2)+(PN+PN)*T(Nl) 

DO  15  1=2, N 

C 1  MOVES  THE  COMPUTATION  ALONG  THE  REMAINING  GRID  POINTS  ON  THE 

C BOTTOM  HORIZONTAL  LINE.  RN=M*T  AT  THESE  POINTS  IS  COMPUTED  HERE. 

C IS  THE  BASIS  FOR  THE  INDICES  OF  Al,  A2,  AND  RN. 

PE=A1(I+1  )+Al( 1  +  2) 

PN=A2( I+N)+A2< I+NN) 

PW=A1 ( I )+Al( 1+1) 

PS=A2( I )+A2( I+N) 

PO=-(PE+PN+PW+PS) 

C FOR  I=N  THE  VALUE  OF  RN(N)  IS  INCORRECT  AND  IS  CORRECTED  BELOW. 

C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

RN( I )=PW*T{ 1-1 )+PO*T< I )+PE*T< 1  +  1 ) ♦ ( PN+PN ) *T ( I+N) 
15  CONTINUE 

C RN=M*T  AT  THE  LAST  GRID  POINT  ON  THE  BOTTOM  HORIZONTAL  LINE  IS 

C CALCULATED  HERE  SEPARATELY  FROM  THE  INTERMEDIATE  POINTS  ON  THE 

C LINE  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  COMPONENT  PE  OF 

C M.  THIS  IS  DONE  BY  SUBTRACTING  PE*T(N+1)  FROM  THE  RN(N)  JUST 

C CALCULATED. 

RN(N)=RN(N)-PE*T(N+1)+PW*T(N-1) 
C N1=N+1,  N21=N**2-2*N+1 

M1=N1 

M2*N 

DO  30  J=N1,N21,N 
r j  MQVFS  THF  COMPUTATION  UP  TO  THE  NEXT  HORIZONTAL  LINE  OF  THE 
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C GRID.  FOR  EACH  J  RN=M*T  IS  CALCULATED  AT  THE  FIRST  GRID  POINT,  THE 

C INTERMEDIATE  POINTS,  AND  THE  LAST  GRID  POINT  OF  THE  GIVEN 

C HORIZONTAL  LINE.  Ml,  M2  AND  J  ARE  THE  BASES  FOR  THE  INDICES  OF  Al 

C A2,  AND  RN.  RN=M*T  FOR  THE  FIRST  GRID  POINT  OF  THE  HORIZONTAL  LINE 

C CALCULATED  HERE  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING 

C COMPONENT  PW  OF  M. 

Ml«Ml+3 

M2=M2+1 

PE»A1(M1)+A1(MH-1) 

PN»A2(M2)+A2(M2+N) 

PW=A1(M1-1)+A1(M1) 

PS=A2(M2-N)+A2(M2) 

PO=-(PE+PN+PW+PS) 
C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

RN(J)=PS*T( J-N)+PO*T(J)+(PE+PE)*T(J+l)+PN*T( J  +  N) 

J1=J+1 

J2=J+N-1 

DO  25  I=J1,J2 

C I  MOVES  THE  COMPUTATION  ALONG  THE  REMAINING  GRID  POINTS  ON  THE 

C GIVEN  HORIZONTAL  LINE  AND  RN=M*T  FOR  THOSE  POINTS  IS  COMPUTED 

C HERE.  Ml,  M2  AND  I  ARE  THE  BASES  FOR  THE  INDICES  OF  Al,  A2  AND  RN 

M1=M1+1 

M2=M2+1 

PE=A1(M1)+A1(M1+1) 

PN=A2(M2)+A2(M2+N) 

PW=A1(M1-1 )+Al (Ml) 

PS=A2(M2-N)+A2(M2) 

PO=-(PE+PN+PW+PS) 

RN( I )=PS*T< I-N)+PW*T( I-l)+PO*T( I )+PE*T( 1+1 )+PN*T( I+N) 
25  CONTINUE 

C THE  COMPUTATION  OF  RN=M*T  AT  THE  LAST  GRID  POINT  ON  THE  GIVEN 

C HORIZONTAL  LINE  IS  DONE  HERE  SEPARATELY  FROM  THE  OTHER  GRID  POINT 

C ON  THE  LINE  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  COMPONENT 

C PE  OF  M.  THIS  IS  DONE  BY  SUBTRACT  I NG  PE*T (J2  +  1 )  FROM  THE  RN(J2) 

C JUST  CALCULATED. 

C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

RN(J2)=RN(J2)-PE*T( J2+1)+PW*T( J2-1  ) 
30  CONTINUE 

L=J2+1 

C L  MOVES  THE  COMPUTATION  UP  TO  THE  FIRST  GRID  POINT  ON  THE  TOP 

C HORIZONTAL  LINE  AND  RN=M*T  AT  THAT  POINT  IS  COMPUTED  HERE.  Ml,  M2 

C AND  L  ARE  THE  BASES  FOR  THE  INDICES  OF  Al,  A2  AND  RN.  NOTE  THE 

C COMPONENT  PN  OF  M  VANISHES  ALONG  THIS  LINE  AND  PW  VANISHES  AT 

C THIS  POINT. 

Ml=Ml+3 

M2=M2+1 

PE  =  A1  (Ml  )+AKMl+l) 

PN=A2(M2)+A2(M2+N) 

PW  =  A1  (Ml-D  +  Al  (Ml) 

PS=A2(M2-N)+A2(M2) 

PO=-(PE+PN+PS+PW) 
C«***#THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

RN(L)=(PS+PS)*T(L-N)+P0*T(L)+(PE+PE)*T(L+1 ) 

L1  =  L  +  1 

DO  45  K=L1,NSQ 

C K  MOVES  THE  COMPUTATION  ALONG  THE  REMAINING  GRID  POINTS  ON  THE 

C TOP  HORIZONTAL  LINE  AND  RN  =  M*T  AT  THOSE  POINTS  IS  COMPUTED  HERE. 
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C Ml,  M2  AND  K  ARE  THE  BASES  FOR  THE  INDICES  OF  Alt  A2  AND  RN. 

M1*M1+1 

M2-M2+1 

PE=A1(M1)+A1(M1>1) 

PN=A2(M2)*A2(M2*N) 

PW«AKMl-l)+Ai(Ml) 

PS*A2(M2-N)+A2(M2) 

PO«-(PE+PN+PW+PS) 
C FOR  K=NSO  THE  VALUE  OF  RN(NSO)  IS  COMPUTED  BELOW. 

IF<K.EQ.NSQ>  GO  TO  40 
C*****JHIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

RN(K)=(PS+PS)*T(K-N)+PW*T(K-1)+P0*T(K)*PE*T(K+1) 

GO  TO  45 

c the  COMPUTATION  RN=M*T  AT  THE  LAST  GRID  POINT  ON  THE  TOP 

C HORIZONTAL  LINE  IS  PERFORMED  HERE  SEPARATELY  FROM  THE  OTHER  GRID 

C POINTS  IN  ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  OF  THE 

C COMPONENTS  PE  AND  PN  OF  M. 

C*#**#THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 
40  RN(K)=(PS+PS)*T(K-N)+(PW+PW)*T<K-1)+P0*T(K) 
45  CONTINUE 

RETURN 

END 
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SUBROUTINE  PRODV( N. S »NSQ . A2 » W, X .NSQP2N ) 
IMPLICIT  REAL*8< A-H.O-Z) 
DIMENSION  S(NSQ).X(NSQ).A2(NSQP2N) 
C 

C CALCULATE  THE  PRODUCT  ( S-( V-W*I ) *X ) .  COMPUTATION  OF  THE 

C COMPONENTS  IS  DONE  IN  DOWN-TO-UP  AND  LEFT-TO-RIGHT  ORDER  TO 

C TAKE  ADVANTAGE  OF  THE  VANISHING  OF  COMPONENTS  OF  V  ALONG  THE 

C BOTTOM  OR  TOP  LINES  OF  THE  GRID. 

c PS.PN.AND  PO  DEFINE  THE  LOWER  DIAGONAL.  UPPER  DIAGONAL.  AND 

C MAIN  DIAGONAL  ( SOUTH .NORTH ,OR IGI N )  OF  V. 

C 

DO  2  I=1,N 

C 1  FOLLOWS  THE  LOWER  HORIZONTAL  LINE.  FOR  EACH  I,  THE  COMPONENTS 

C OF  THE  PRODUCT  ARE  COMPUTED  ALONG  THE  CORRESPONDING  VERTICAL  LINE 

c THE  FIRST  COMPUTATION  FOR  EACH  I,  AND  THE  LAST.  ARE  SPECIAL 

C CASES  OCCURRING  WHEN  PS  AND  PN  ARE  ZERO. 

C i    MOVES  THE  COMPUTATIONS  FROM  LEFT  TO  RIGHT  ALONG  THE 

C (BOTTOM)  HORIZONTAL  LINE. 

IPN=I+N 

PS=A2( I )+A2( IPN) 

PN=A2( IPN)+A2( IPN+N) 

PO=-(PS+PN) 
C SAVE  X(I)  TO  USE  AGAIN 

X1=X( I ) 
C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

X( I  )=S( I  )-(  (PO-W)*X< I )+(PN+PN)*X(IPN)  ) 

C THE  VARIABLE  OF  THIS  LOOP  STEPS  UP  EACH  VERTICAL  LINE  IN 

C ORDER  FROM  LEFT  TO  RIGHT.  THE  LOOP  VARIABLE  IS  THE  INDEX  OF  X. 

C HENCE  IPN  IS  THE  INDEX  OF  X  AT  THE  LOWEST  POINT  AND  LI  IS  THE 

C INDEX  OF  X  AT  THE  HIGHEST  POINT. 

L1=NSQ-N-N+I 

DO  1  J=IPN.L1»N 

PS  =  PN 
C JPN  IS  THE  INDEX  OF  A2 

JPN=J+N 

PN=A2(JPN)+A2( JPN+N) 

PO=-(PS+PN) 

X2=X( J) 

X(J  )=S( J)-(PS*Xl+(PO-W)*X( J)  +  PN*X(JPN) ) 

X1  =  X2 

1  CONTINUE 

C THIS  SEGMENT  COMPUTES  THE  COMPONENTS  ALONG  THE  UPPER 

C HORIZONTAL  LINE. 

PS  =  PN 

C K2  IS  THE  X  INDEX 

C L1  =  N**2-2*N  +  K 

K2=L1+N 
C K3  IS  THE  A  INDEX 

K3=K2+N 

PN=A2(K3)+A2(K3+N) 

PO=-(PS+PN) 
C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

X(K2)=S(K2)-( (PS  +  PS)*Xl  +  (PO-W)*X(K2)  ) 

2  CONTINUE 
RETURN 
END 
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SUBROUTINE  HH1 ( N, Al , NSQMN, W,H,NSQ, NSOP2N ) 

IMPLICIT  REAL*8(A-H,0-Z> 

DIMENSION  H(NSQ),A1(NSQP2N) 
C 

C THIS  SUBROUTINE  COMPUTES  THE  SOLUTION  OF  (H+W*I)X«S  BY  GAUSSIAN 

C ELIMINATION,  THE  MATRIX  (H+W*I)  IS  TRIDIAGONAL  WITH  THE  MAIN 

C DIAGONAL  AND  THE  TWO  ADJACENT  DIAGONALS  ON  EITHER  SIDE  OF 

C THE  MAIN  DIAGONAL. 

C COMPUTATION  IS  DONE  IN  THE  FOLLOWING  MANNER:  FIRST  FORWARD 

C ELIMINATION  IS  PERFORMED  TRANSFORMING  THE  TRIDIAGONAL  MATRIX 

C INTO  AN  UPPER  DIAGONAL  MATRIX  (SUBROUTINE  HH1).  THEN  (SUBROUTINE 

C PP1)  THE  RIGHT  HAND  SIDE  VECTOR  IS  TRANSFORMEDt  THEN  THE 

C BACKWARD  SUBSTITUTION  IS  PERFORMED  TO  SOLVE  THE  RESULTING 

C UPPER  TRIANGULAR  SYSTEM.  H  IS  SET  EQUAL  TO  ZERO  ALONG  THE 

C RIGHTMOST  VERTICAL  LINE  OF  THE  GRID,  WHICH  RESULTS  FROM  THE 

C FACT  THAT  THE  RIGHT  ADJACENT  DIAGONAL  VANISHES  AT  THESE  POINTS. 

C 

DO  1  K=N,NSQ,N 

H(K)=0.0 

1  CONTINUE 

C NEXT  THE  COMPUTATION  IS  PERFORMED  ON  THE  FIRST  VERTICAL  LINE  OF 

C THE  GRID  SINCE  H  AT  THESE  POINTS  DEPENDS  ONLY  ON  THE  DIAGONAL 

C ELEMENT  AND  THE  RIGHT  ADJACENT  DIAGONAL  ELEMENT,  AND  DOES  NOT 

C DEPEND  ON  ANY  PREVIOUS  H. 

C PW,PE,AND  PO  DEFINE  THE  LEFT  DIAGONAL,  RIGHT  DIAGONAL,  AND 

C MAIN  DIAGONAL  ELEMENTS  OF  THE  MATRIX  H. 

NTOP=NSQMN+l 

M=l 

DO  2  I=l,NTOP,N 

C 1  IS  THE  H  INDEX  AND  MOVES  THE  COMPUTATION  UP  THE  FIRST  VERTICAL 

C LINE.  IM  IS  THE  Al  INDEX. 

IM=I+M 

PW=A1 ( IM-1 )+Al( IM) 

PE=A1 ( IM)+A1( IM+1) 

PO=-(PE+PW) 
C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

H(I )=-(PE+PE)/(PO+W) 

M  =  M  +  2 

2  CONTINUE 

C THE  COMPUTATION  OF  H  I S  NOW  PERFORMED  IN  LEFT-TO-RIGHT  AND  DOWN- 

C TO-UP  ORDER;  I.E.  H  IS  CALCULATED  ALONG  THE  REMAINING  GRID  POINTS 

C ON  THE  FIRST  HORIZONTAL  LINE  STARTING  WITH  THE  FIRST  AND  I  MOVES 

C THE  COMPUTATION  UP  TO  THE  NEXT  HORIZONTAL  LINE. 

N1=N+1 

M=l 

I  1=2 

INM1=N-1 

DO  5  1=1, N 

PW  =  A1  (  I+M)+A1  (  I+M+l  ) 

DO  4  J=l I t INM1 
C J  IS  THE  H  INDEX  AND  Jl  IS  THE  Al  INDEX 

J1=J+1 

PE  =  A1  (Jl  )+Al ( Jl  +  1 ) 

f>0»-(  pf-f  pw) 

H(J )=-PE/ ( PQ  +  W+PW*H( J-l  )  ) 

^--^ 
<*    CONTINUE 
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II=II+N 

INM1=INM1+N 

M=M+N1 

CONTINUE 

RETURN 

END 
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SUBROUTINE  PP1 ( N, Al ,NSQMN, W,H ,NSQ,X,NSQP2N) 
IMPLICIT  REAL*8(A-H,0-Z) 
DIMENSION  H(NS0),X(NSQ),A1(NSQP2N) 

C 

C this  SUBROUTINE  CONTINUES  THE  GAUSSIAN  ELIMINATION  PROCESS 

C BY  CALCULATING  THE  RECURSION  FORMULA  FOR  THE  RIGHT  HAND  SIDES; 

C I.E.  IT  CALCULATES  P(M)=X(M)  FOR  M«1,...,N  BY  THE  FOLLOWING 

C FORMULA: 

C 

C P(M)=X(M)=(X(M) -PW*X ( M-l ) ) / ( PO+W+PW*H ( M-l ) ) 

C WHERE 

C P(l)=X<l)=X(l)/(PO+W) 

c 

C AND  WHERE  PW,  PE,  AND  PO  ARE  THE  LEFT  DIAGONAL,  RIGHT  DIAGONAL, 

C AND  MAIN  DIAGONAL  ELEMENTS  OF  H  AT  POINT  M. 

C COMPUTATION  IS  DONE  IN  A  LEFT-TO-RIGHT  AND  DOWN-TO-UP  ORDER 

C OVER  THE  GRID.  FOR  EACH  HORIZONTAL  LINE  P(M)=X(M)  IS  CALCULATED 

C AT  THE  FIRST  POINT  AND  THEN  AT  THE  REMAINING  POINTS  IN  ORDER  TO 

C TAKE  ADVANTAGE  OF  THE  FACT  THAT  P(M)  DOES  NOT  DEPEND  ON  A  PREVIOUS 

C H  AT  THE  FIRST  POINT  AND  DOES  DEPEND  ON  A  PREVIOUS  H  AT  THE 

C REMAINING  POINTS.  THE  I  LOOP  MOVES  THE  COMPUTATION  UP  TO 

c THE  NEXT  HORIZONTAL  LINE  OF  THE  GRID. 

C 

INN  =  1 
11=2 
IN  =  N 

M=l 

DO  2  1=1, N 

C COMPUTATION  FOR  THE  FIRST  POINT  ON  THE  GIVEN  HORIZONTAL  LINE 

C IS  PERFORMED  HERE.  INN  IS  THE  X  INDEX  AND  IM  IS  THE  Al  INDEX. 

IM=INN+M 

PW  =  A1  (IM-1  )+AK  IM) 

PE=A1( IM)+A1 ( IM+1) 

PO=-(PE+PW) 

X( INN)=X( INN)/(PO+W) 

DO  1  K=II,IN 

C THIS  LOOP  PERFORMS  THE  COMPUTATION  AT  THE  REMAINING  GRID  POINTS 

C ALONG  THE  GIVEN  HORIZONTAL  LINE.  K  IS  THE  X  INDEX  AND  KM  IS  THE 

C Al  INDEX. 

KM=K+M 

PW  =  PE 

PE=A1 (KM)+A1(KM+1 ) 

PO=-(PE+PW) 

KK=K-1 
C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

IF(K.EO.IN)GO  TO  4 

X(K  )=(X(K)-PW*X(KK)  )/(PO+W+PW*H(KK)  ) 

GO  TO  1 

4  X(K)=(X(K)-(PW+PW)*X(KK) ) / ( PO+W+ ( P W+PW ) *H ( KK ) ) 

1  CONTINUE 
INN=INN+N 
J  I  =  1  I+N 
IN= IN  +  N 

M  =  M+2 

2  CONTINUE 

C       THIS  SECTION  COMPUTFS  THE  BACKWARD  SUBSTITUTION  FOR  THE  GAUSSIAN 
IMINATION  PROCEDURE;  I.E.  IT  COMPUTES: 


c  89 

C X(N)=P(N)=X(N) 

C X(M)*P(M)+H(M)*X(M+1)*X(M)+H<M)*X(M+1)  FOR  M=N-1,...,1 

C 

C XI  SAVES  X(NSQ)  TO  USE  AGAIN 

X1=X(NSQ) 

I1»NSQ-1 

DO  3  1=1,11 
C KK  IS  THE  DECREASING  INDEX  FOR  X  AND  H 

KK=NSQ-I 

X(KK)=X(KK)+H(KK)*X1 
C XI  SAVES  X(KK)  TO  USE  AGAIN 

X1=X(KK) 
3  CONTINUE 

RETURN 

END 
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SUBROUTINE  PRODH< N, S,NSQ, Al ,NSQMN, W, X, NSQP2N ) 

IMPLICIT  REAL*8( A-HfO-Z> 

DIMENSION  S(NSQ),X(NSQ>,A1(NSQP2N) 
C 

C CALCULATE  THE  PRODUCT  ( S- (H-W*I ) *X ) •  COMPUTATION  OF  THE  COMPONENT 

C IS  DONE  IN  LEFT-TO-RIGHT  AND  DOWN-TO-UP  ORDER  TO  TAKE  ADVANTAGE 

C OF  THE  VANISHING  OF  COMPONENTS  OF  H  ALONG  THE  LEFTMOST  AND 

C RIGHTMOST  VERTICAL  LINES  OF  THE  GRID. 

C PW,PE»AND  PO  DEFINE  THE  LEFT  DIAGONAL,  RIGHT  DIAGONAL*  AND  MAIN 

C DIAGONAL  < WEST , EAST, OR IGEN )  OF  H. 

C 

M  =  l 

NTOP=NSQMN+l 

DO  2  l=l,NTOP,N 

C 1  FOLLOWS  THE  LEFTMOST  VERTICAL  LINE  OF  THE  GRID.  FOR  EACH  I  THE 

C COMPONENTS  OF  THE  PRODUCT  ARE  COMPUTED  ALONG  THE  CORRESPONDING 

C HORIZONTAL  LINE. 

C THE  FIRST  COMPUTATION  FOR  EACH  I,  AND  THE  LAST  ARE  SPECIAL 

C CASES  OCCURRING  WHEN  PW  AND  PE  ARE  ZERO.  THIS  FIRST  SECTION 

C PERFORMS  THE  COMPUTATION  AT  THE  FIRST  POINT  ON  THE  GIVEN 

C HORIZONTAL  LINE. 

c i     iS  THE  x  INDEX  AND  IM  IS  THE  Al  INDEX.  XI  SAVES  X(I)  TO  USE 

C AGAIN 

IM=I+M 

X1=X(I ) 

PW  =  A1  (IM-D+AH  IM) 

PE=A1( IM)+A1( IM+1) 

PO=-(PE+PW) 
C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

X(I )=S(I )-( (PO-W)*X( I )+(PE+PE)*X(I+l ) ) 

C THIS  SECTION  PERFORMS  THE  COMPUTATION  ALONG  THE  INTERMEDIATE 

C POINTS  ON  THE  GIVEN  HORIZONTAL  LINE. 

11=1+1 

I2=I+N-2 

DO  1  J=I1,I2 

C J  IS  THE  X  INDEX  AND  JM  IS  THE  Al  INDEX.  X2  SAVES  X(J)  TO  USE 

C AGAIN. 

JM=J+M 

X2  =  X(J  ) 

PW  =  PE 

PE=A1( JM)+A1( JM+1 ) 

PO=-(PE+PW) 

X(J  )=S(J)-(PW*X1+(P0-W)*X( J)+PE*X(J+1 ) ) 

X1  =  X2 

1  CONTINUE 

C THIS  SECTION  PERFORMS  THE  COMPUTATION  AT  THE  LAST  POINT  OF  THE 

C GIVEN  HORIZONTAL  LINE. 

L= I  2  +  1 

LM=L+M 
C L  IS  THE  X  INDEX  AND  LM  IS  THE  Al  INDEX. 

PW  =  PE 

PE=A1 (LMJ+A1 (LM+1 ) 

PO=-(PE+PW) 

^♦♦♦♦♦this  change  is  required  for  the  neumann  problem. 
/'l)=s(l)-(( pw+pw)*x1+(p0-w)*x(l) ) 

m  =  m  +  ;> 

2  CHNTINUF- 
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RETURN 
END 
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SUBROUTINE  HH2 < N, A2 , NSOMN , W ,H , NSQ » NSQP2N ) 

IMPLICIT  REAL*8( A-H»0-Z> 

DIMENSION  H(NSQMN)»A2(NSQP2N) 
C 

C THIS  SUBROUTINE  SOLVES  (V+W*I)X«S  BY  GAUSSIAN 

C ELIMINATION.  SINCE  V  IS  A  TRI-DUGONAL  MATRIX  CONSISTING  OF  A 

C MAIN  DIAGONAL  AND  TWO  OFF  DIAGONALS  AT  A  DISTANCE  N  FROM  THE 

C MAIN  DIAGONALt  GAUSSIAN  ELIMINATION  CAN  BE  PERFORMED  DIRECTLY 

C ON  THE  MATRIX, 

C FIRST  THE  FORWARD  ELIMINATION  PROCESS  IS  PERFORMED 

C TRANSFORMING  THE  TRI-DIAGONAL  MATRIX  INTO  AN  UPPER  DIAGONAL  MATRIX 

C (SUBROUTINE  HH2 )  THEN  (SUBROUTINE  PP2)  THE  RIGHT  HAND  SIDE  VECTOR 

c IS  TRANSFORMED  IN  ACCORDANCE  WITH  THE  UPPER  DIAGONAL  MATRIX;  I.E. 

C THE  P(I)=X(I)  ARE  CALCULATED.  THEN  THE  BACKWARD  SUBSTITUTION 

C (SUBROUTINE  PP1  )  IS  PERFORMED  TO  GET  THE  RESULT. 

C ALONG  THE  FIRST  HORIZONTAL  LINE  OF  THE  GRID  THE  VALUES  OF 

C THE  UPPER  DIAGONAL  COMPONENTf  H,  DO  NOT  DEPEND  ON  ANY  PREVIOUS 

C VALUES  OF  H.  PS»PN  AND  PO  DEFINE  THE  LOWER  DIAGONALt  UPPER 

C DIAGONAL  AND  MAIN  DIAGONAL  ELEMENTS  OF  V. 

C THE  COMPUTATION  IS  DONE  IN  A  DOWN-TO-UP  AND  LEFT-TO-RIGHT 

C ORDER  TO  TAKE  ADVANTAGE  OF  THE  VANISHING  OF  COMPONENTS  OF  V  ALONG 

C THE  BOTTOM  AND  TOP  LINES  OF  THE  GRID. 

C 

DO  2  1=1, N 

C i  FOLLOWS  THE  LOWER  HORIZONTAL  LINE.  FOR  EACH  I  THE  VALUES  OF  H 

C FOR  THE  TRANSFORMED  UPPER  DIAGONAL  MATRIX  ARE  COMPUTED  ALONG  THE 

C CORRESPONDING  VERTICAL  LINE.  I  MOVES  THE  COMPUTATIONS  FROM  LEFT  TO 

C RIGHT  ALONG  THE  (BOTTOM)  HORIZONTAL  LINE.  I  IS  THE  H  INDEX  AND 

C IPN  IS  THE  A2  INDEX. 

IPN=I+N 

PS=A2(I )+A2( IPN) 

PN=A2(IPN)+A2( IPN+N) 

PO=-(PS+PN) 
C*****THIS  CHANGE  IS  REOUIRED  FOR  THE  NEUMANN  PROBLEM. 

H(I )=-(PN+PN)/(PO+W) 

C J  STEPS  UP  EACH  VERTICAL  LINE  IN  ORDER  FROM  LEFT  TO  RIGHT.  J  IS 

C THE  H  INDEX  AND  JPN  IS  THE  A2  INDEX. 

L1=NSQ-N-N+I 

DO  1  J=IPN,LltN 

PS  =  PN 

JPN=J+N 

PN=A2( JPN)+A2( JPN+N) 

PO=-(PS+PN) 

H( J )=-PN/(PO+W+PS*H( J-N) ) 

1  CONTINUE 

2  CONTINUE 

C THE  VALUE  OF  H  FOR  THE  TRANSFORMED  UPPER  DIAGONAL  MATRIX  IS  ZERO 

C ALONG  THE  TOP  HORIZONTAL  LINE  OF  THE  GRID  SO  IT  NEED  NOT  BE 

C COMPUTED  AND  STORED. 

RETURN 

END 
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SUBROUTINE  PP2 ( N, A2 , NSQMN, W»H ,NSQ » X , NSQP2N ) 

IMPLICIT  REAL*8( A-H»0-Z) 

DIMENSION  H( NSQMN) , X ( NSQ > * A2 < NSQP2N ) 
C 

c THIS  SUBROUTINE  CONTINUES  THE  GAUSSIAN  ELIMINATION  BY  CALCULATING 

c THE  RECURSION  FORMULA  FOR  THE  RIGHT  HAND  SIDES;  I.E.  IT  CALCULATE 

c THE  RIGHT  HAND  SIDES  P(M)*X(M)  BY  THE  FOLLOWING  FORMULA: 

c PIM)xX(M)  FOR  M*lt...»N 

C P(M)»X(M)=<X(M)*PS*X(M-N))/(PO-PS*H<M-N) )  FOR  M«N+1,..,N 

c WHERE  PStPN  AND  PO  DEFINE  THE  LOWER  DIAGONALf  UPPER  DIAGONAL 

C AND  MAIN  DIAGONAL  ELEMENTS  OF  V  AT  POINT  M. 

c COMPUTATION  IS  DONE  IN  A  DOWN-TO-UP  AND  LEFT-TO-RIGHT  ORDER 

c OVER  THE  GRID.  FOR  EACH  VERTICAL  LINE  P(M)=X(M>  IS  CALCULATED  AT 

c THE  FIRST  POINT  AND  THEN  THE  REMAINING  POINTS  IN  ORDER  TO  TAKE 

C ADVANTAGE  OF  THE  FACT  THAT  P(M)  DOES  NOT  DEPEND  ON  A  PREVIOUS  H 

C AT  THE  FIRST  POINT  OF  ANY  VERTICAL  LINE. 

C 

DO  2  1=1, N 

C 1  FOLLOWS  THE  LOWER  HORIZONTAL  LINE.  FOR  EACH  I  THE  VALUES  OF 

c p-x  ARE  COMPUTED  ALONG  THE  CORRESPONDING  VERTICAL  LINE.  I  MOVES 

C THE  COMPUTATIONS  FROM  LEFT  TO  RIGHT  ALONG  THE  (BOTTOM)  HORIZONTAL 

C LINE.  I  IS  THE  INDEX  OF  P  (P=X)  AND  IPN  IS  THE  A2  INDEX. 

IPN=I+N 

PS=A2( I )+A2( IPN) 

PN=A2( IPN)+A2( IPN+N) 

PO=-(PS+PN) 

X(I )=X( I )/(PO+W) 

C J  STEPS  UP  EACH  VERTICAL  LINE  IN  ORDER  FROM  LEFT  TO  RIGHT.  J  IS 

C THE  X  INDEX  AND  JPN  IS  THE  INDEX  OF  A2. 

L1=NSQ-N+I 

DO  1  J=IPN»L1»N 

PS  =  PN 

JPN=J+N 

PN=A2(JPN)+A2( JPN+N) 

PO=-(PS+PN) 
C*****THIS  CHANGE  IS  REQUIRED  FOR  THE  NEUMANN  PROBLEM. 

IFIJ.E0.L1 )GO  TO  4 

X(J)=(X( J)-PS*X( J-N) )/(PO+W+PS*H( J-N) ) 

GO  TO  1 
4  X<J)=(X( J)-(PS+PS)*X< J-N) )/(PO+W+(PS+PS)*H( J-N) ) 

1  CONTINUE 

2  CONTINUE 

C THIS  SECTION  OF  THE  SUBROUTINE  PERFORMS  THE  BACKWARD  SUBSTITUTION 

C FOR  THE  GAUSSIAN  ELIMINATION  PROCESS;  I.E.  IT  COMPUTES: 

C 

C X(M)=P(M)=X(M)  FOR  M=N, . . . t N**2-N+l 

C X(M)=P(M)+H(M)*X(M+N)=X(M)+H(M)*X(M+N)  FOR  M=N**2~N t . . . ♦ 1 

C 

I 1=NSQMN+1 

DO  3  I=l»NSOMN 
C KK  IS  THE  DECREASING  INDEX  FOR  X  AND  H. 

KK=I1-I 

X(KK)=X(KK)+H(KK)*X(KK+N) 

3  CONTINUE 
RETURN 
END 
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0.991878550  01 
_Q._829C61700_01 
0.683186540  01 
0.551332210  01 
0.43089739D  01 
0.320930910    01 


0.538329710 
0.479161910 
0.414559960 
0. 34530928 D 
0.275364320 
0.20676433D 


0.476564090-04 
0. 303064  930-04 
0.171179260-04 
0.833717690-05 
C. 336751150-05 
0.  107327410-05 
0.253160600-06 
0.40647076D-07 
0.397306370-C8 
0.202569950^09 
0. 42^619900-11 
0.260637310-13 


0.223247960  01 
_0.  141764760    01 

0.800733660  00 
.j0_._389.99 18  50    QO 


00 
00 
00 
00 
00 
iO_Q_ 


0.157523590  00 
0.502C50570-01 


0.  14417918D  00 
0.915887010-01 
0.517322  03  0-01 
0__2519  5827D-0i 
0.101769730-01 
0.32435495  0-02 


0.118421970-01 
■O.JL  901371  5  0-02 
0.185848350-03 
0.947602110-05 
0.200914090-06 
0.122441560-08 


0.765077260-03 
0.  122840070-03 
0.  12006922D-04 
0.61220829D-06 
0.  1293  02  870-07 
0.79141226D-10 


C. 227595720-15 
0  .292358730-15 


0.645156150-11 
0.632494060-1  I 


0. 505463010-13 
0.166317100-13 
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ITERATION       0 


ITERATION  I 
_LT£RAIi.QN 2_ 

ITERATION  3 
_IJLERAimN__4_ 

ITERATION   5 
ITERATION   6 


MAXIMUM  NORMALIZED  RESIDUAL 
0.666666670-01 


ITERATION       7 

JLERAJLIONl 8_ 

ITERATION      9 

J_I£RAI.IQN_JLQ_ 
11 
_12_ 


ITERATION 
J_I£BATJ_ON_ 


ITERATION 
JJEKATION 
ITERATION 
ITERATION 
ITERATION 


13 


15 

_1_6 

L7 

JLIE_RAT_iU.M_lJi. 


ITERATION 
I  TE  RATION. 
ITERATION 
_IXE  RATION 
I TERATION 
ITERATION, 


19 
_20_ 


21 
2_2_. 
21 
24 


ITERATION  25 

-LTERJAT_IUI\L_26. 

ITERATION  27 
_LI£RAJLJON_28 

ITERATION  29 
-1_LEKAI_IUN_3.QL 


ITERATION    31 
JJERALI_QN_J2_ 


0.  774839750- 
J3_f_1294Q36  20- 

0.A09652390- 
JL.2135_4547D- 

0.268492610- 

0.484439610- 


L2    NCKM 
0-0 


RATIO 


03 
0  4_ 

06 
06 
07 

08 


0.227C7515D 
-Q_.J0_9l32l50j 
0.33  1200160- 
0_._1 1 39  96  17J> 
0.101195940- 


0.722432800-02 
JD.6C799540D-03 
0.12  82  71540-03 
0.498937000-04 
0.936444  720-05 
0.332202750-05 


■08 
08 
09 
Q9_ 
10 
0.21839503D-11 


0.12C555250 
JLa_£_Q  l.C  8  3  63  0  ■ 


0.102915060- 
0.264337650 
0.355649  80  0- 
0.117  20  26  10- 


0.717648800- 
-0  _.3_4  0  43  2_3  3  D- 
0.913163020- 
0.389883400- 
0.107633010- 
^.JJ^_9287_8_0r 
101907850- 
_1116A06  2D- 
1096  36810- 
_1Q4  341040- 


0 
_0 

0 
_0_ 

0.807247560- 
JL._S1A5_9_431I>: 


11 

12 

12 

13_ 

14 

i_4_ 
15 
15 
16 
16_ 
16 
I6_ 


0.52734296D-06 

JLl2_58151_310-06 
0. 805263900-07 

-P.  •  JL1 03  3  5  2  6  0  -  0  7_ 
0.504483430-08 
0.171819690-08 
0.  302390000-09" 

_0.  132522740-09 
0.254808720-10 

_0_.„7619  5  591D-11 
0.27822733D-11 
0.136527380-11 


0.212719510-12 
_0_.  95897955  0-13 
0.355086330-13 
0.21699  89  30-13 
0.180295930-13 
_0.  178457760-13 


0. 112499390- 
_C_._78_720_9  5  00- 


16 
l_6 
16 
16 
17 

17_ 

16 

17 


0.174681090-13 
_0_._i7193l310-13 
0.171073140-13 
0_.  16446736  D- 13 
0.164683530-13 
0.164271250-13 
0.16  250  3730-13 
0_._158  757030-1 3 
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PAX  IMiJM    NHRMAI  T/FD    RESIDUAL         LI    VFCTUR    NORM         L2    NORM    RATIO 

I  0.797638270-03  0.132027410    02      0.79557158D   00 

— 2 0*21 9 C 98 8 8D_r.O  3 0, 16331 76  3D_0 1 Q.  6 76  7 7 0 38 D   00_ 

3  0.140652440-03  0.564218360    01      0.578310120    00 

— 4 0  _X2.06.661 50-0  3__* 0.45Q3108  90    01  _  0.  4922  7654  0_QlO_ 

5  0.10169786D-03  0.359253310    01      0.405728460    00 

— 6 0. 841913240- OA Q-__e.8ji___.L3  2  0  01      0.330909970   pp 

7  0.677108250-04  0.230288240   01      0.26614257D    00 

8 _U3133  79  890-04       0.  181339620„  Ol_„  0*2  10604830    Q0_ 

9                      0.33077995D-04                         0.  11 1 14725D    01      0.130136180    00 
_ 1__ 0_.  201  132A20j^0A Q._66786606_0_  00 0^786700.9  7  D-iOl 

11                      0.119375860-04                         0.391223090    00      0.463211190-01 
-12 -L^__-_i_66  320-Q5 -_-£____A_6840    00      0. 2  6447227D-0 ] 

13                      0.381611070-05  0.121712130    00      0.145383  110-01 

-14 _U______£l^l_5J_rJ_5 0.6 3930894LH0  1  _0.  766663370-02 

15                     0.103131300-05  0.320013560-01      0.385200880-02. 

^6 •CljlA.S33.52_62  0ji.06. 0,15159.0400-01      0.  1831241  90-02_ 

17                      0.12128716D-06                        0.353128240-02      0.436659050-03 
-^ Q.?75.9_33A.6J0_r__J 0.789096120-03      0.  9691  5161  D-04 

19                      0. 576003310-08  0.159898130-03      0.197599490-04 

-2& 0.108395390-08 Q_.  294  3  893.1 0-04  „  0.  365741  11  D-05 

21                      0.182199790-09  0.481213610-05       0.600641480-06 

-22 Cl*_26_69_422  3^-L0 . C__691  i_3433D--g6_  _0x866l53860-07_ 

23                      0.335122460-11  0.853194880-07      0.107299830-07 

-2A Q_-.3517  5  573JUZLL2 0_y_3_8268148L)-08      0.111347660-08 

25                      0.28986685D-13  0.72  3561600-09      0.908414150-10 

-2& 0.22C26265D-14 O.A973  36520-10  _0.  578305550-1 1 

27  0.149838090-14  0.162947840-10      6.2  74531 710- I 2~ 

28 0^1.15762730-  .14. 0 .JA3JI 2  3  5  0- 1 0_      0.  120402  32D-l_3_ 

29                      0.103671200-14  0. 1 53913900-10      0.135199540- 13 

-^ Q.1338^650Jl_J-4 0.168981130-10      0.  18434904D- 1  3 

31  0.153786330-14  0.169494620-10       0.171655260-13 

32  0.142934340-14  __0.  160742210- 10      0.133213100-1 3__ 
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ITERATION 

0 

MAXIMUM  NORMALIZED  RESIDUAL 
0.675675680-01 

L2  NORM  RATIO 
0.0 

ITERATION 
ITERATION 
ITERATION 
I  TERATION 

I 
? 
3 

4 
5 
6 

0. 943201420-02 
0.21729C41D-02 
0.43241693D-03 
0.88390354D-C4 

0.58071917D  00 
0.384511640  00 
0.33300406D  00 
0.28885845D  OC 

ITERATION 
ITERATION 

C.81952496D-04 
0.76944366D-04 

0.272475510  00 
0.25716231D  00 

ITERATION 
ITERATION 

7 

8 

9 

10 

11 

12 

0.65054589D-04 
0.57914335D-04 

0.195799160  00 
0.149466790  00 

ITERATION 
ITERATION 

0.51634272D-C4 
-  C. 408757680-04 

0.13478101D  00 
0.  12170359D  00 

ITERATION 
ITERATION 

C. 416565630-04 

0.37906C480-04 

0.116314710  00 
0.111237850  OC 

ITERATION 
ITERATION 

13 
14 
15 
16 
17 
18 

C. 327552990-04 
0.276181640-04 

0.915899870-01 
0.75465564D-01 

ITERATION 
ITERATION 

0.261786300-04 
C.22081712D-04 

0.698832840-01 
0. 647777970-01 

ITERATION 
ITERATION 

0.226677680-04 
0.211568260-04 

0.624851140-01 
0.60303519D-01 

ITERATION 

ITERATION 

19 

20 

C. 154642440-04 
C.14249C730-C4 

0.413547370-01 
0.2860C5790-01 

ITERATION 
ITERATION 

21 
22 
2  3 

24 

0.114200C6D-04 
0.871719260-05 

0.24833732D-01 
0.216134940-01 

I TERATION 

ITFRATION 

0. 791038830-05 
0.674848950-05 

0.2038 402 6D-0 I 
0.19245750D-01 

ITERATION 
ITERAT  ION 
ITERATION 
ITERATION 

25 

26 
27 
23 

29 
30 

0.552276400-05 
0.479452C30-05 
0.41538816D-C5 
0.325878950-C5 
0.328332930-05 
0.296451740-05 

0.147042170-01  . 
0.112552510-01 
0.101484940-01 
0.916555210-02 

ITERA  T  ION 
ITERATION 

0.87539842D-02 
0.837677770-02 

ITERATION 
ITERATION 

31 
32 

0. 252968720-05 
0.212497900-05 

0.689941040-02 
0.563688060-02 
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: MAXIMUM    NCRMALIZEQ    RESIDUAL 

I  0.494798920-02 

-2 Q_,  2843356  1  0-02 

3  0.195144450-02 

.4 0_a  16  62  2  63  6_Dr_0_2_j 

5  0.14742072D-02 

-A 0.133137820-07 


LI    VECTOR    NORM         12    NORM    RATIO 


9 
JLQ_ 


11 
12_ 


13 
_L4_ 


15 
_L6_ 

17 
JLfl. 


0.  120948990-02 
._£)__.  .1 098  544  2  0-  0  2 

0.995361680-03 
_  JO  ,.905  090820-03 

0.111484520-02 

C.139C733Q0-Q2 


C.20C499770-02 
0.2818046 10-0 2 
0.389338490-02 

XL.  5 l27  2606  50-0  2 
0.7017  05990-02 

XL^_9 136301  50-02 


19 

_2.0_ 

21 

22 


0.113329830-01 
Q._15C3  72_90O-0l 
0.  188353290-01 
,0.232626260-01 


23 

_2_4_ 


0.28382900D-01 
0.342741270-01 


25 
26 


27 

23 
29 


0.410170270-01 
0.486  7463  00_-_0_l_ 
0. 572637160-01 
0.  66  75266  70-31 
0.769726740-01 
0.876125280-01 


0.758910990 

_0_.  57833  31.60 

0.471624960 

0.39746038D 

0.34C83210D 

0.295145260 

0.258891210   02 

.0.22  7617060    02 

0.202  593350    02 

_0_._18_400  377  0    02 

0.1727116  20    02 

0.  171604630  02 


02       0.95094974D    00 

.02 0_.  920531390  00 

02  0.89465124D  00 
02  0.8 700  3 394  0  00 
02  0. 845064 C70  00 
02       0. 818631020    00 


0.  790042  760  00 
0. 758386110  00 
0.72295505D  00 
0.682953190  00_ 
0.63  7549090  00 
0.535354460    00 


0.18139187D  02 

JL.202918170  02 

0.236487040  02 

_0.  2 81 2  568 10  02 

0.339072860  02 

0.413396430  02 


0.504574450  02 
_0.61240J340  02 
0.735960170  02 
0_.JT75_672_16O  02 
6.103362  010  0  3 
0.  121357690    03 


0.527115090  00 
0.461026670  00 
0.383303870  00 
0.3U32799D  00 
0.234642720  00 
0.165285670  00 
0.11338167D  00 
0.906478230-01 
0. 95895432 D-01 
0.113267370  00 
0.  134655730  00 
0.  15953446D    00 


0.142056750  03 

JD.  165944500  03 

0.193189530  03 

_0. 223804360  03 

0.2  569  86350  0  3 

0.290985500  03 


0.  188666050 
0.222440490 
0.260939320 
0. 303912070 


31 
32 


0.98 13  13470-0  1 
0.107725950    00 


0.322851940    03 
0.348697680    03 


00 
00 
00 
00 
00 
00 

0.445377730   00 
0.485431830    00 


0.350490380 
0.398764720 
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MAXIMUM    NORMALIZED    RESIDUAL         L2    NCRN    RATIO 

ITERATION      0 0 . 67  1  76877Q-C  1 Q.Q 

ITERATION       1  C.e92877130-C2  0.580917300    00~ 

-IHRATjCN 2 C.15A21A550-C2 0.38  75189<5D    00 

ITERATION      3  C  .25728335D-C3  ~07T359T5TlD~ 6d" 

_1_TERATI0N      A 0.  830A70  15D-0A 0.2917921  CD    OC 

ITERATION      5  0 .820  273700-CA  "  T72T53l5'35D~0C" 

ITERATION      6 C  .7776  10AQD-CA 0.26C0AQ700    00 

ITERATION       7  0.659789600^  0TT9827108C    00 

ITERATION       8 0.574627970-C4 0.15166685D    CO 

ITERATION       9  C.  A9  6C6  8  6  6D-CA  ~07T369Tr2C'D~00" 

JLI E R AJUO N_ jm C.A1A92 CJB6 D-CA  0,1  21725. 8 6D    00 

ITERATION     11  C.A161857?D-CA  "T0riTS3C5C90~00" 

ITERATION     12 0_.378A627CD-CA  0.113183130    OC 

IT  ERA  TICK    n  0 .  3 32836580-CA  0.9330e597~f>oT~ 

ITERATION    1A  Q. 2  7 9C956 C  D^C A 0.77C02A71D-01 

ITERATION    15  C.26C13916D-CA  ~CT71 36  3~1 520-0  V 

JJER  ATICN    16 C226C  8  A5  6D-C  A g  .661897580-01 

ITERATION    17  C  .  2288AC520- CA  ~0~  63  37 18  C  6  0  HJT~ 

ITERATION     18 C  .  2  10C2  1  79D-CA 0  .6  1660336D-Q1 

ITERATION    19  C. 1 58C325CD-CA  0. A 2356229D-C I 

JTERATICN    20  C  .  1A070  1  78D-CA 0  .  293A5259D-C1 

ITERATION    21  0  .  1082306  1D-0A  07255Tsl9T6^0T" 

ITERATICN    22 0  .8A6672860-C5 0.2223  1700D-Q1 

ITERATION    23  C  .  78  1 50AC30-C5  ~T720V8_1387C^0T~ 

ITERATION    2A C  .685  182  79C-C5  0.  198  18A96D-C  1 

ITERATION    25  C  .  5668C A 360-C5  0.  1  516AC86C-C1 

JJER ATJT  N_  2b 0  .  A8398599D-C5_  _  0.116277A50-01 

ITERATION    27  C  .A05201  37C-C5  "  ^Tl  CVsTC  SW- CT~ 

U  ERA_11C_N_2_8 C  .33255  9^  2D-  C  5 0  .9A8666060-C? 

ITERATION    29  0 .33381 8A2D-C5  ~^0T9C7b~A6ClC-02~ 

ITERATION    30 0  .  298893  33D-C5 0.867776750-02 

ITERATION    31  0.2628559C0^C5  0. 7156A31AD-02 

ITERATION    32 C  .2  18A8  1  7  1D-C5  0  .59C76955D-02 
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5 
6_ 


—MA.X1HUM    NORMALIZE: J    RESIDUAL 
1  0.564897290-02 

_2 0_*  29  5  4 .4  3_6 5  0-  0 .2_ 

3  0.225026750-02 

0  , 1_9£JL_Q  5000- C2_  ; 

0.  180U133D-02 
0.16  5742  7JLQr£L2 


LI    VECTOR     NORM  L2    NORM    RAT|n 


0. 693666390 
Q ,  5 327  79  7a  D 
0.434 A 8336  0 
_0j  3669 700 ID 
0. 3 1 58 31 78 D 
0.274832070 


9 
.JO. 
11 
_L2_ 


0.  1526  30330-02 
_C-._i3t970702D-02_ 
0.126531500-C2 
0.115776280-02 
0. 12  33  77210-02 
-Q...13123  082  0-02 


0.240801760 
_0_,_212504  370 

0.189578620 
0.17539019D 
0.170659230 


13 

15 
_16_ 

17 
_L8_ 


0  .179791520 


0.185037900-02 
_  _Q_.  25  77  3939  0-0  2_ 

0.350730460-02 
_J). . 4 6  0 3 V.  18  6 0-0 2 

0.6197  30880-02 

0.808013300- 


19  0.104750090-01 

_  2  0 _0_.  1319  930  70-0 1 

21  0.163947000-01 

JL2^„  _0.210  166770-01 

23  0.258303940-01 

_2it OJ.  14  1133  50-01 

25  0. 37334343 J-OL 

26 0.451510100-01 

27  C.  533607850-01 

2d  0^.623  735920-0  1 

29  0.719957750-01 

_1Q O.Bia345750-01 


0.204563040  02 
0.244791330  02 
0.302629430  02 
_0.3„802  063  50_02 
0.477167450  02 
0.593006450    0? 


0.728576570 
0.381267470 
0.104851150 
.0.122891130 
0.142281030 
0.163263320 


0.186038400 
0.211098400 
0.238659410 
0.268506630 
0.299893300 
0.331496220 


02  0.95195615D  00 
0 ?___ D j. 9  21972520  00 
02      6.896417590   00 

02 0.87208958D_00_ 

02      0.84  74  11140    00 

.Qi 0.821354290  00 

02  0.79311204D  00 
02_0,761960  39D  00 
02   0.727186630  00 

0_2 _0_.  6 88046970^00. 

02  0.643  74633D  00 
02  0. 593458430  00 
0.53643826D  00 
0._4_L23Q4710_00 
0.40152653  0  00 
0.326007650  00 
0.249570580  00 
0.17828852D  00 
02   0.121090990  00 

02  0.9014 5 140  0-01 

03  0.90285282D-0r 

03 0.106_98389jD  J)0_ 

03  0.128591640  00 
03  0.153133350  00 
03  0.  181304910  00 
03  0.213611390  00_ 
03  0.2  5010  8560  00 
03_  0^2903  9578  0  00 
03  0.333408030  00 
03       0.377115530    00 


0.913102260-01 
0,996320100-01 


0.361179630   03 
0.386320870    03 


0.418308430    00 
0.4527218  70    00 
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MAXIMUM    NORMALIZED    RESIDUAL         L2    NORM    RATIO 
_TTERATION      0  0.900000000    00  0.0 

I TERATION       I  0.36565756D    CO  0.9996  71210    0~0~ 

ITERATION       2  0. 8489971 7D-CU  0.99887916D    00 

ITERATION      3  0.165282650-01  I)7997Tr6~3  8  D~0"0" 

ITERATION       4  0. 1 77350 75D-01 0.99634342D    00 

ITERATION      5  0  .37  3  1 72  I A I'D-  02  0T9957  C3~5~2D~ 00" 

ITERATION       6 0  .309260560-02 0.  9950782  9D    00 

ITERATION      7  0. 809555210-02  0.992611170    00 

ITERATION       8  0  .  933506220-C2 0.990061UD    00 

ITERATION      9  6.242126280-02  "0Y9890 375lD~00" 

ITEPATION     10  CU1 18377  890- C2  0.9880 249 7D    OC 
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Results  of  ADI  Algorithm  on  Heterogeneous  Model  with  Random  No.    Sub- 
regions  with  Neumann  Boundary  Conditions 
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